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PART I INTRODUCTION 


In the past decade there has been much activity in the development of computational methods 
for the analysis of unsteady transonic aerodynamics about airfoils and wings. Advances have 
paralleled developments in steady computational fluid dynamics (CFD) with a lag of 
approximately five years 1 due to the additional requirement of time-accuracy. Also 
contributing to this time lag is the sheer number of calculations required to perform flutter 
analyses, a primary application of unsteady CFD. Figure 1 , taken from the specification 


document for U.S. military aircraft.2 illustrates significant features which must be addressed 
in the treatment of computational unsteady transonic aerodynamics. On the plot of equivalent 
airspeed versus Mach number, lines of constant altitude are straight lines through the origin 
with decreasing altitudes represented by lines with steeper slopes. An airplane's flight 
envelope is typically set by the maximum limit speed and a typical flutter boundary curve, 
characterized by the flutter airspeed gradually dropping to a minimum in the transonic speed 
range followed by a rapid upward rise. The ability to predict this minimum, termed the 
transonic flutter dip, is of great importance in design, since the flutter boundary must be shown 
by a combination of analysis and test to be outside the flight envelope by a margin of at least 15 
percent in equivalent airspeed, i.e. the flutter boundary must be outside the dashed line 
boundary in fig. 1. Subsonic linear unsteady aerodynamic theories have been quite successful in 
predicting this flutter boundary for Mach numbers up to 0. 6-0.7 but linear theory does not 
account for the effect of aerodynamic shape or maneuvering conditions upon unsteady airloads at 
transonic speeds. At these Mach numbers linear analysis has been used with more or less 
success depending upon the severity of local transonic effects. The occurrence of flutter within 
the flight envelope usually leads to structural failure and loss of the vehicle, highlighting the 
necessity of careful validation of computational methods intended for use in this area. This is a 
key difference in the utilization of steady and unsteady computational methods which should be 
clearly understood. 


Transonic Flow Phenomena 

It will be helpful to distinguish the main features of steady transonic flow in order to 
organize the discussion of unsteady aerodynamics. Figure 2, from ref. 3, indicates various 
regions of transonic flow development for the NLR 7301 airfoil, a 16-percent thick cambered 
supercritical-type section. With increasing Mach number and moderate angle-of-attack, the 
upper surface becomes critical between M = 0.4-0. 7 with the first shock forming at an 
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increase of approximately 0.1 in Mach number. Pearcy et al. 4 have classified several types of 
flow separation which may occur. For conventional airfoils the typical pattern, termed type A, 
involves the growth of a local separation bubble induced by boundary layer separation at the 
shock foot, spreading rapidly to the trailing edge as Mach number increases. This condition is 

often accompanied by unsteady phenomena such as buffet and aileron buzz 3 . The steep aft 
pressure gradients of modern airfoils, such as the NLR 7301, can lead to an alternate pattern, 
termed type B, in which separation progresses from the trailing edge towards the shock. Figure 
2 illustrates this type B separation, with fully separated flow aft of the shock occuring along the 
line of maximum lift. Note the small "shock free" design condition occuring over a small 
isolated range of lift coefficient and Mach number just prior to the onset of trailing edge 

separation. Tijdeman 3 notes the flow conditons in the region between the onset of trailing edge 
separation and fully separated flow are very sensitive to Reynolds number and the location of 
transition from laminar to turbulent flow. 

Figure 3 shows a similar diagram, derived from ref. 5. of attached, mixed, and separated 
flow regions for a complete aircraft at free stream Mach numbers between 0 and 2.0. In region 
I, the flow is predominantly attached. To obtain optimum performance and to avoid the* diag 
penalty associated with flow separation, design cruise conditions for aircraft typically are 
located in region I, near the boundary of region II (mixed flow). 

As speed and/or angle of attack increase, a transition region of mixed flow (region II of fig. 

3) is encountered. For rigid structures, this region is typified by the onset of localized regions 
of flow separation which may exhibit significant aerodynamic unsteadiness. For realistic 
flexible structures, the aeroelastic response of the structure interacts with the airflow to 
induce much more complicated situations. For instance, structural vibrations can cause the 
flow to alternately separate and reattach at flow conditions where a rigid structure would 
support attached flow. The associated highly unsteady aerodynamic loading can interact with the 
structure to cause unusual aeroelastic phenomena which may restrict the vehicle flight 
envelope. 

With further speed and/or angle of attack increases which may be encountered under 
maneuvering conditions, stable separated flow conditions emerge (region III of fig. 3). Leading- 
edge vortex flows and shock-induced vortex flows are of this nature. At still higher angles, 
vortex bursting in the vicinity of the aircraft can cause severe buffeting. Within such regions 
the flow is highly unsteady and accurate computations will require careful attention to 
turbulence modeling. - 
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While predictive methods for attached flows are reasonably well developed, the picket fence 
in fig. 3 emphasizes the difficulty in predicting aeroelastic phenomena in the mixed and 
separated flow regions. It also symoblizes novel features that are being encountered in 
transonic flutter testing. Modem high performance aircraft are capable of maneuvering at 
transonic speeds, leading to a much enlarged parameter space that must be considered in flutter 
analysis and resting. Wing/store ioading, fuseiage interference, angle-of-attack, wing shape 
and wing sweep all must be considered, and the traditional flutter boundary parameterization of 
dynamic pressure at flutter versus Mach number may need to be augmented to adequately 
describe aeroelastic stability boundaries. For instance, flutter tests give some indication that 
these additional parameters affect the detailed aeroelastic stability condition near the flutter 
boundary. Thus, the pickets of the fence in fig. 3 represent possible regions of low damping or 
instability that might be encountered. 

Farmer et al. 6 provided early test results documenting the effect of airfoil shape upon 
flutter boundaries. Figure 4 shows their comparison of flutter boundaries for two structurally 
and geometrically similar wings of the same planform. The supercritical wing was a reduced 
stiffness model of the modified TF-8A wing while the conventional wing had a symmetrical 
section. The two wings had leading-edge sweep angles of 44.5 degrees. Design cruise Mach 
number was 0.90 for the conventional wing and 0.99 for the supercritical wing. The 
supercritical wing was shown to have a 25 percent lower minimum flutter dynamic pressure 
near Mach 1.0 where type II mixed flow would be expected. Current transonic computational 
methods are beginning to address this important area which will be a key topic for 
computational aeroelasticity in the future. Other reports of aeroelastic model tests relevant to 

this area are; single mode flutter of a low aspect ratio wing studied by Erickson 7 , supercritical 
wing flutter tests performed at the NLR 8 . 9 and torsional buzz of aeroelastic wings tested at the 

FSAE10. 

Figures 5-8 illustrate four types of aeroelastic response which have been encountered and 
which offer challenges for computational methods. The four cases illustrate problem areas 
encountered near the boundaries of aircraft flight envelopes, as operating conditions change 
from high speed, low angle conditions to lower speed, higher angle conditions. The nonclassicai 
aeroelastic response observed on the DAST ARW-2 wing model 11 , fig. 5, is a region of high 
dynamic response at nearly constant Mach number which was encountered at dynamic pressures 
well below those for which flutter was predicted. The motion is of the limit-amplitude type and 
the response is believed to be associated with flow separation and reattachment over the 
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supercritical wing (type II flow). 

Figure 6 illustrates wing/store limited amplitude oscillations experienced by modern, high 
performance aircraft under various loading and maneuvering conditions at transonic Mach 
numbers. Such oscillations can result in limitations on vehicle performance. The conditions 
for which this response occurs appear to be near the onset of type II mixed flow. The response 
typically increases for maneuvering flight conditions. 

Dynamic vortex-structure interactions causing wing oscillations have been observed on a 

bomber type aircraft 12 for high wing sweep conditions during wind-up turn maneuvers, fig. 7. 
The flow involves the interaction of the wing vortex system with the first wing bending mode 
and occurs over a wide Mach number range (0.6-0.95) at angles of attack of 7-9 degrees. 

At higher angles, interaction of forebody and wing vortex systems with aft vehicle componets 
results in vortex-induced buffet loads, illustrated in fig. 8. The figure shows the operating 
conditions for which tail buffet may occur on a high performance fighter. Buffet of horizontal 
tails can occur at intermediate angles of attack and is a result of the vortex system encountering 
the horizontal tail lifting surface. As angle of attack increases, the location of vortex bursting 
moves upstream in the wake. Loss of lift is associated with the burst location reaching the 
vicinity of the aircraft, and vertical tail surfaces located in such regions can experience severe 
dynamic loads. 


Historical Perspective 

This field received an initial impetus in the mid-1970's from three sources: Tijdeman's 3 
pioneering experimental work on transonic unsteady pressure measurements, Magnus and 
Yoshihara's demonstration of key transonic flow features for an airfoil with an oscillating 
flap 13 and the introduction of an economical transonic finite-difference solution algorithm by 

Ballhaus and Goorjian 14 . Ballhaus 15 gives a survey of the field from this period. The AGARD 
Structures and Materials Panel Subcommittee on Aeroelasticity has selected experimental 
unsteady pressure data sets and defined two- and three-dimensional Standard Aeroelastic 

Configurations 16 - 17 to provide reference computational test cases for the development and 
validation of improved computational methods. The data sets were obtained from rigid models 
undergoing pitch and control surface oscillations and includes both conventional and 

supercritical airfoil geometries 1 8,19.20. | n addition to these data sets, Sandford et al.2i 

summarizes a series of unsteady pressure tests made at NASA Langley and Tijdeman 22 presents 
a much used data set for a fighter wing configuration. 
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Computational methods have been pursued at a number of differing levels of physical 
approximation to the flow equations. Magnus and Yoshihara^ 3,23,24 used an explicit algorithm 
to solve the Euler (EE) equations. Steger and Baiiey25 reported a significant early application 
to the problem of aileron buzz using an implicit approximate factorization solution algorithm 
for the navier-Stokes (NS) equations. Chyu and his coauthors26.27 have pursued further 
applications of derivatives of this code. Most of the nonlinear unsteady computations to date 
have been made by solving the potential equations, both with and without interacted viscous 
effects. For example, the alternating-direction implicit (ADI) algorithm embodied in the 
LTRAN2 code of Ballhaus and Goorjianl 4 enabled efficient solutions of the two-dimensional flow 
frequency transonic small disturbance (TSD) potential equation through the use of large time 
steps. Extensions of this ADI algorithm have been widely used by many researchers. A semi- 
implicit form of the ADI algorithm is used in the 3-D XTRAN3S code28,29 developed for the 
aeroelastic analysis of wings. Other TSD and Full Potential (FP) equation codes are described in 
refs. 30-35. There is a growing trend, especially for steady flows, towards use of the Euler 

equations rather than the potential equations. Euler equation codes treating 2-D oscillating 

% 

airfoils are reported in refs. 36-40 while SaImond4l and Belk42 show results from 3-D Euler 
codes. 

Over this same time period, several experimental investigations of periodic aerodynamic 
flows about rigid airfoils have been reported. McDevitt 4 3,44 documented such conditons for a 
very narrow range of Mach number of an 18 percent thick circular arc airfoil and Levy45 

reproduced the effect with calculations from a NS code. Subsequently, Mabey 46 studied these 
oscillations for circular arc airfoils with thicknesses of 10-20 percent. References 47 and 48 
give details for a 14 percent circular arc airfoil. Related information regarding the interaction 
of unsteady airloads caused by transitional boundary layers with structural oscillations is given 
by Mabey et al 49 . Another class of separation-induced periodic flow problems, vortex shedding 
about rigid cylinders and airfoils at high angle-of-attack, has been studied using NS codes for a 
variety of Reynolds numbers in refs. 50-52. 

Unsteady aerodynamics has been the theme of four recent AGARD conferences 5 3-56 whose 
proceedings contain a wealth of information. Survey papers focusing upon computational 
requirements and resources are given by Peterson 6 ^ and McCroskey et al. 66 . Summary papers 
of the 1984 and 1985 AGARD conferences are given by Mykytow59 and by Mabey and 
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Chambers 60 . The latter reference makes recommendations regarding computational and 
experimental methods for unsteady flow phenomena and draws particular attention to the need to 
pay careful attention to the nature of shock motions. The periodic oscillations about circular 
arc airfoils are recommended as benchmark computational cases for all time-dependent 
transonic viscous flow theories. Zwaan 61 surveys aeroelastic problems in transonic flow while 

Deiwert 62 reviews the numerical simulation of unsteady interactive flows. Finally, Mabey 63 
gives a review of pertinent experimental research on time-dependent aerodynamics. 

Experimental Data Sets 

In this section, the airfoil geometries and wing planforms which have been most frequently 
studied are summarized. In addition to the AGARD standard configurations, several other model 
tests have been popular for comparison with computational results. Figures 9 and 1 0 show the 

profiles and planforms of the 2-D 1 6 and 3-D 1 7 AGARD configurations, respectively. Data sets 
for all of these configurations except the 6 percent parabolic arc, DO A1 and MBB-A3 airfoils 
are given in refs. 18, 19, and 64. Tables 1 and 2, from ref. 1, tabulate selected references for 
these and other configurations in which comparisons of experimental and calculated unsteady 
pressures are given. The entries are grouped by the equation level of the physical modeling used 
for the calculations. The references are not exhaustive but are an attempt to indicate 
publication of significant experimental/computational comparisons or new capability. 

The first three airfoils in Table 1 are conventional airfoils with 6, 10, and 12 percent 

thickness ratios. Tijdeman 3 tested the NACA 64A006 airfoil with an oscillating quarter-chord 
trailing-edge control surface. Interpretations of these tests 3 have provided insights into the 
underlying mechanisms of unsteady transonic flows. Tijdeman indentified three types of shock 
motion, denoted type A, B, and C. In type A shock motion, the shock wave remains distinct 
during the oscillation cycle, with a periodic variation of shock location and shock strength. In 
type B shock motion, the shock wave weakens and disappears during a portion of the cycle, 
generally during the forward propagation of the shock along the surface. For type C motion, the 
shock wave on the airfoil remains distinct and propagates forward along the airfoil chord and off 
the airfoil leading-edge. 

Davis and Malcolm 66 tested the NACA 64A010A airfoil for pitching oscillations. Two cases 
from this test have been widely studied: a case with a moderate shock wave at M = 0.8 and a = 0 
degrees and a case with steady shock-induced separation at M - 0.8 and a = 4 degrees. The NACA 
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0012 airfoil, tested by Landonis, differs from the other entries in Table 1 in that it was tested 
for larger dynamic pitching amplitudes and for transient ramping motions making it suitable 
for dynamic stall computational studies. McDevitt and Okuno 66 have reported measurements of 
periodic shock-induced oscillations for this airfoil. 

Data sets for the 16 percent thick supercritical NLR 7301 airfoil are given by both 


Tijdeman and Davisi 8 and the shock-free condition for this supercritical airfoil has been a 
challenging computational case. The 8.3 percent thick MBB A-3 asrfos! has been tested by 


Zimmerman 87 and represents a less severe supercritical airfoil computational case. Other 
supercritical airfoils tested for oscillatory motions or exhibiting unsteady behavior are: a 12 
percent thick airfoil tested for pitching, heaving and flap rotation by den Boer and Houwink 88 , 
the RA16SC1 airfoil tested by ONERA 69 , and the cryogenic test of a supercritical SC(2)-0714 


airfoil by Hess et al 70 . Reference 68 reported large dynamic responses of airloads on the 
supercritical airfoil for both oscillating and static motions at type II flow conditions and 
introduced the concept of "aerodynamic resonance." Similar periodic shock-induced oscillations 


are reported for the RA16SC1 airfoil 89 . 

Tests of rigid circular arc airfoils have been reported by McDevitt et al. 43 , McDevitt 44 , 
Mabey 46 and Mabey et al. 47 . References 43 and 44 give details of tests of an 18 percent thick 
airfoil for Reynolds numbers of 1 million to 17 million, covering laminar to fully developed 
turbulent flows. The wind tunnel walls were contoured to approximate the inviscid stream- 
lines over the airfoil at M = 0.775. Periodic unsteady airflows were observed over a narrow 
Mach range whose extent depended upon whether Mach number was increasing or decreasing. 

For increasing Mach numbers, oscillations occurred for 0.76 < M < 0.78 while for decreasing 
Mach number the range was wider, 0.73 < M < 0.78. The frequency of the oscillations was 188 
±3 Hz (reduced frequency k = 0.48 based upon semi-chord). Mabey 48 studied similar periodic 
flows for a series of circular arc airfoils ranging in thickness from 10 to 20 percent at 
Reynolds numbers of 0.4-0.6 million. In ref. 47, further investigations on a larger 14 
percent thick biconvex wing at Reynolds numbers of 1-7 million is reported. Two necessary 
criteria evident from the experimental results for the existence of the periodic unsteady flow 
are given: thickness/chord ratio greater than 12 percent and local Mach number upstream of 
the terminal shock wave in the range 

1 .24 < M < 1 .40 


McDevitt 44 identifies the predominant shock motion for the 18 percent thick airfoil as type C 
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whereas Mabey et al. 47 argue that it is type B motion. 

The smaller number of entries in Table 2 reflects the situation regarding 3-D testing in 
that there are fewer experimental data sets widely available and fewer comparisons of 
experimental and calculated results have been published. Tijdeman22 tested a model of the F-5 
fighter wing including external tanks and stores. This wing has an aspect ratio of 2.98, a taper 
ratio of 0.31 and a leading edge sweep of 32 degrees. The relatively thin wing section, a 
modified NACA 64A004.8. has made this a popular computational case since it is well within the 
capability of TSD codes. Transonic and low supersonic test conditions are available. Of the 
AGARD Standard Configuration models shown in fig. 6, the NORA model is the most extensively 
tested. It is a model of the Mirage F-1 horizontal tail which has been tested in four European 

wind tunnels 17 - 18 . 

The AGARD rectangular wing and the RAE Wing A model have symmetric airfoil 
sections 1 7 - 18 -20 whereas the ZKP wing and LANN wing have supercritical airfoil sections 17 - 19 . 
Additional models tested for oscillatory pitching are the NASA Rectangular Supercritical Wing 
(RSW) model 71 - 72 and the RAE AGARD tailplane model 64 . The former had a 12 percent 
supercritical airfoil section while the latter had a NACA 64A010A section, the same as one of 
the AGARD 2-D configurations. 

Also included in Table 2 are references to several other published comparisons with 
experimental data. These cases are of interest since the models were aeroelastic and some 
comparisons of experimental and computed transonic flutter boundaries (or aeroelastic 
response) are given. Isogai gives comparisons for a high aspect ratio supercritical transport 

wing in ref. 34 and for the supercritical wing flutter model of Farmer et al. 6 in ref. 74. 

Bennett et al. 75 give static aeroelastic comparisons for an aspect ratio 10.3 supercritical wing 
which was extensively instrumented for unsteady pressure measurements 76 . Finally, 
Guruswamy and Goorjian 77 present calculations for a rectangular parabolic arc flutter model. 

Computational Methods 

A variety of fluid dynamic flow models is available to address unsteady aerodynamic 
computations. The choice of an appropriate method calls for assessment of the difficulty of the 
aerodynamic problem being addressed. Type I flows, fig. 3, include one of the most important 
aeroelastic analysis conditions, cruise at high dynamic pressure. Classical linear aeroelastic 
analysis has been primarily focused upon this condition. The transition from type- 1 to type II 
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conditions can occur due to aircraft maneuvering with little decrease in dynamic pressure. 
Thus, aeroelastic response and stability of aircraft operating in type II flows can be quite 
important although they have only recently been brought within the range of computational 
methods. 

Computational methods available for unsteady aerodynamic computation include; the 
classical (linear) small disturbance poieniial equaiion (CSD), nonlinear potential equation 
(both Transonic Small Disturbance, TSD, and Full Potential, FP, Euler equations (EE) and 
Navier-Stokes equations (NS). 

Issues which have been central to unsteady CFD have been the choice of implicit versus 
explicit algorithms, the stability of alternative solution algorithms and the treatment of 
computational grids. Explicit schemes are simple to code and easily vectorizable but are limited 
in allowable time step by the stability limit imposed by the signal propagation time over the 
smallest grid cell. Faced with the requirement of maintaining time-accuracy throughout the 
entire field for aeroelastic computations, this easily leads to excessive computation times, 
especially for viscous flow calculations where a very fine mesh near the surface is required to 
^resolve the boundary layer. The alternative implicit solution algorithms thus are favored and 
attention must be given to their relative stability and accuracy characteristics. Grid generation 
for unsteady problems in which the body boundary moves, such as for an oscillating control 
surface or an aeroelastic deformation, raises new issues over those involved in steady flows. To 
maintain accuracy, the body-conforming grid must be realigned with the body at each time step. 
Schemes for accomplishing this have been studied as well as the necessity of moving the grid at 
all. When body motions are small with perturbations mainly normal to the surface, imposing 
boundary conditions on the mean surface location may be an acceptable approximation. Finally, 
the nature of unsteady calculations means that the solution is not allowed to achieve a steady- 
state and thus the dynamic response of numerical calculations on the computational grid is more 
important. For example, grid cell stretching in the near and far field will affect the 
computational impedance of the grid for unsteady calculations. 

In the following sections the physical flow models will be described in order of decreasing 
complexity, to be followed by discussions of typical results illustrating progress for the 
various flow modeling levels. Since the historical trend is for computational methods to mature 
most rapidly for the simpler flow models, this will lead us naturally from classical linear 
results back to the Navier-Stokes equations. Along the way, capabilities such as the ability to 
treat geometric complexity will be discussed. 
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PART II FLUID DYNAMIC FLOW MODELS 
Navier-Stokes Equations 

Rumsey and Anderson78 give the thin-layer approximation to the Reynolds averaged Navier- 
Stokes equations for two-dimensional flow. In the thin-layer approximation viscous terms are 
resolved in a layer near the body where viscous terms in the direction along the body, are 
neglected and only terms in ti, normal to the body are retained. The equations are written in 
generalized coordinates and conservation form; 
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The curvilinear generalized coordinates (^, r\) correspond to the coordinates parallel and 
normal to the body surface, respectively and are related to Cartesian coordinates (x, y) via the 


transformation 
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£ = (x, y, t), rj = (x, y, t), x = t ( 4 ) 

Note that the transformation is time-dependent, allowing the grid to move to follow body motion 
and giving rise to grid metric terms such as t|;p in eq. (3). The contravariant velocities along 
the and £ coordinate directions are 


(5) 

V = r| x u + ri y v + r| t 

while the pressure is 

P = (Y-D[e-^p(u 2 + v 2 )] (6) 

The state vector Q represents the density, momentum and total energy per unit volume. The 
Jacobean of the transformation is J, defined as; 


a(x, y) 


(7) 


The equations are nondimensionalized by the freestream density p~ and soundspeed a«. The shear 
stress and heat flux terms are defined in tensor notation as; 


M 


0 u. d u. 


X x x =-^ UK-^ + ^ + X 
x = x - Re, d x. d x. 

L j i 


» j 


d u, 

, 8..] 
d x, *j 


•_ _ , M . ,3<» 2 ) 

L J 

1 Re L Pr(y-l) d x i 


( 8 ) 

(9) 


Re L = ' 


p q L 



M 


( 10 ) 



14 


In (3), b is defined as 

x i 


b = UT 


* q 

*i x j 


( 11 ) 


Stokes hypothesis for bulk viscosity, X. + 2^/3 = 0, and Sutherland's law for molecular 
viscosity, 

p. = p>|T = (T/T ) 3/2 [(T +c)/(T + c)] (12) 

oo oo oo 


are used, with T~ = freestream temperature = 460° R, and c = Sutherland's constant = 198.6° 
R. 

Boundary conditions are applied explicitly. No slip, adiabatic wall conditions, as v/ell as 
zero normal pressure gradient conditions are applied on the body 

u = v = 0 (13a) 


3 p _ d (a 2 ) 
di\ 3ri 


(13b) 


In the farfield, a quasi-one-dimensional characteristic analysis is used to determine boundary 
data. For turbulent calculations, turbulence modeling such as the algebraic eddy viscosity model 

of Baldwin and Lomax 79 is required. 


Euler Equations 

For sufficiently large Reynolds numbers the major effect of viscosity is confined to a thin 
viscous boundary layer near the surface of a solid body. As a consequence, the inviscid portion 
of the flowfield may be solved independently of the boundary layer. The reduced set of equations, 
termed the Euler equations, are obtained by dropping both the viscous terms and the heat 
transfer terms from the Navier-Stokes equations. Reference 80 contains details and discusses 

consequences of these assumptions. Anderson et al.® 1 present the three-dimensional Euler 
equations in generalized coordinates for a moving grid mesh in a form analogous to eqs. (1-7). 
With the generalized coordinates ^ and £ prescribed parallel and rj normal to the body surface 
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the equations are 
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The pressure p and the contravariant velocities U, V, and W are given by straightforward 
extensions of eqs. (5) and (6). The grid transformation metric variables are given by; 
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In ref. 81 the boundary conditions are applied explicitly. On the body, the contravariant 
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velocity normal to the body is set to zero 

V = V/|gradr|| =0 (18) 

and pressure, tangenial velocities and density are determined by extrapolation from the 
interior. In the farfield, a quasi-one-dimensional characteristic analysis is used to determine 
boundary data. 

While the Euler equations do not treat viscous and heat transfer effects, entropy and 
vorticity effects allow treatment of flows with strong shock waves (which generate entropy) 
and moderate-to-high angle vortex dominated flows. Even though viscosity is eliminated, no 
explicit Kutta condition enforcing smooth flow from the trailing edges of lifting surfaces is 
required since "numerical" viscosity generated by finite difference solutions provides this 
effect. 


Potential Equation 

If it is assumed that the flow is irrotational then the velocity field, V, can be shown to be the 
gradient of a scalar field variable, the velocity potential <t> (see ref. 80 for details) 


V = V<D (19) 

The conservation form of the continuity equation for two-dimensional flow becomes (ref. 82) 

P t + (P^ x ) x + (P^A = 0 (20) 

For barotropic, isentropic flow the momentum and energy equations yield the compressible 
Bernoulli equation from which the density p is determined. 

l 

P = [l+^W-2® VV] 7 '' (21) 

The spatial coordinates, x and z, are normalized by airfoil chord c, and time t is normalized by 
dicJc. Density and <t> are normalized by p~ and a~c. 

Again, solutions are obtained for body fitted generalized coordinates which allow for body and 
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grid motion. Note that here and in the following small disturbance potential equation, the 
coordinates x and z are the freestream direction and the direction normal to the lifting surface 
respectively. The transformation to a body-fitted coordinate system is given by 

£ = £ (x, z, t), C = C ( x * z » 0. x ~ 1 ( 2 2 ) 


The strong conservation form of eq. (20) is maintained by writing the continuity equation as 

<^ +(f 7V (£ rV° <23> 


Equation (21) transforms to 

l 

p = ( 1 + 111 [M* - 2<D x - (U + ^ t )<E>„ - (W + } 71 

Where the contravariant velocities are 


u -s. + <£ + #V (V, + 

(25) 

w -c. + <V,+WV ( £ + $®. 


and 

Bernoulli equation as 


The isentropic pressure coefficient is derived from the compressible 


C p = — — y (p Y * 1) (26) 

yM 

oo 

The airfoil boundary condition that the flow be tangent to the airfoil is satisfied by requiring W 
= 0 on the airfoil. For lifting flows, the shed vorticity is represented by a jump in potential 
across the wake line. For isentropic flows the condition is 


T +<W>r =0 (27) 

T C 

where r is the jump in potential across the wake, <t> u - <l> 1 , and <W> is the average of W above 
and below the wake. In the farfield, the flow is set to freestream conditions 
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O = M x , p = 1 


(28) 


Transonic Small Disturbance Potential Equation 
The Transonic Small Disturbance (TSD) Potential equation is derived from the inviscid 
Euler equations assuming that the flow is a small perturbation of a steady uniform flow, U~, in 

the x direction (see, for instance, ref. 83). The TSD velocity potential function, <>, describes 
the perturbed velocity components u, v, w. 


u 


3<J> 

a? 




(29) 


where the total velocity in the x direction is U~ + u. References 84 and 85 give the modified 
TSD potential equation in conservation form as 


3f o 9fj 9f 2 9f 3 

9t + dx + 9y + Bz 


= 0 


(30) 


where 

f o = - A V B 0x 
f 1= E(l) x + F^ + G(t)y 

(31 ) 

f 2 = + H<t> x <f) y 

f 3=* Z 

Time, t, and the Cartesian coordinates x, y, and z are nondimensionalized by the freestream 
velocity and wing reference chord. 

The coefficients A, B and E are defined as 

A = M 2 , B = 2M 2 , E = 1 - M 2 


(32) 
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Several choices are available for the coefficients F, G and H depending upon the assumptions used 
in deriving the TSD equation29. Briefly, the coefficients are referred to as "NASA Ames" 
coefficients when defined as 

F = -^(Y+1)M 2 , G = 1(y-3)M 2 , H = -(y-1)M 2 (33) 

and are referred io as the “NLR“ coefficienis when defined as 

F = - i [3 - (2 - y)M 2 ]M 2 , G = -Im 2 ,H = -M 2 (34) 

The "classical" TSD coefficients are given by 

F = -^(Y+1)M 2 , G = 0, H = 0 {35) 

and finally the coefficients for the linear potential equation, valid for subsonic and supersonic 
small perturbation flows, are 


F = G = H = 0 (36) 

The TSD equation (29) is distinguished from the higher equation level flow models in that, 
within the small disturbance assumption, the computational grid is not required to move with 

the body since boundary conditions are imposed at the mean plane, usually z = 0 ± . The wing 
flow tangency condition boundary condition is 

= (37) 

where f± (x,y,t) = 0 describes the upper and lower body surfaces. The trailing wake boundary 
conditions are 


[<\+<tg=o 


(38) 


tog = o 


(39) 



20 


where [ - ] indicates the jump in the indicated quantity across the wake. Equation (37) 
enforces the convection of vorticity downstream from the trailing edge and eq. (38) requires 
continuity of the 2 component of velocity, w, across the wake. The pressure coefficient may be 
computed using either linear or nonlinear forms of the Bernoulli equation. The exact nonlinear 
equation is given by eq. (30) where the appropriate density equation for eq. (30) is 

l 

p = [l -^M^(2<{) x + (l)^ + (J)y + <!)^ + 2<j) t )] Y ' 1 (40) 


Alternatively, the linearized Bernoulli relation gives 


c p = - 2( <j ) x + <|>t) (41) 

While eq. (41) is the proper choice based upon formal order of magnitude reasoning, the higher 
order terms in eq. (40) are sometimes not negligible for cases of interest. 

Farfield Boundary Conditions : Two forms of farfield conditions have been used for the 
unsteady TSD equation. Table 3 lists these as "reflecting" and "nonreflecting" boundary 
conditions. These terms are descriptive only and indicate the relative effects of the two sets of 

conditions upon unsteady calculations. The nonreflecting conditions are derived 87 from a 
characteristic variable analysis of solutions to eq. (30) in the farfield. They are implemented 
as first order plane wave conditions and are intended to prevent the reflection of a major 
portion of signals incident upon the boundaries back towards the vicinity of lifting surfaces. 

Note that in the steady state, the time derivatives in eqs. (40 b-g) vanish resulting in simple 
Neumann "reflecting" boundary conditions. Proper grid design is very important for unsteady 
calculations and involves consideration of grid extent, grid stretching and boundary conditions. 

If any one of these factors is not properly treated, spurious unsteady results may be observed. 

It is noted that an alternative farfield boundary condition has been shown to be superior for 
two-dimensional steady calculations. Imposing conditions appropriate to a point vortex of the 

appropriate strength to match lift 87 has yielded solutions of the Euler equations which were 
insensitive to computational grid extent, leading to efficient solutions for small grid extents. No 
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corresponding development for unsteady flows is available. 


Table 3. Farfield boundary conditions for the Unsteady TSD potential equation 



D nflftntrn« 
■ tv/iigvm 

N o n reflect in i g 


Upstream 

<() = 0 

<t> = 0 

(40a) 

Downstream 

<l> x + <l> t = o 

1 -B D v 

2 ( ~C + -7= )( >t + t = 0 
V ^ 

M0b> 
% ~ $ 

Above 

O 

II 

•O-* 

D 

J*t + +z m0 

(40c) 

Below 

II 

o 


( 4 0 d ) 

Right spanwise 

o 

II 


(40e) 

Left spanwise ^ —q 

(for full-span modeling) 


(4 0 f ) 

Symmetry plane 6-0 

Y y 

(for half-span modeling) 

o 

il 

(40g) 


C = E + 2(J> x , D - yj 4 A + B 2 /c 


Linear Small Disturbance Potential Equation 
When the coefficients of eq. (36) are used in eq. (30) the classical linear potential equation 
results. In dimensional form as given in ref. 88 it is 

2 1 2M ? 

VV-Ld) - — <t> - M 2 <b =0 

Y 2 Vtt a Vxt y xx 

a « 


( 43 ) 
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Equation (43) holds for small disturbances from freestream conditions for subsonic and 
supersonic flight. 

Special cases of eq. (43) have been extensively studied and closed form solutions are 
available for some cases. In other cases, asymptotic methods provide insight into the form 
required of numerical solutions. Since these linear solution methods are well calibrated with 
regard to important issues, such as wing flutter, for speeds up to high subsonic Mach numbers, 
it behooves practitioners of unsteady CFD to verify their methods, wherever appropriate, 
against these linear methods. 

The important case of steady state, simple harmonic motion has been most extensively 
studied. With body motion and potential assumed given by 

f* (x, y, t) = f ± (x, y)e ,<ot (44a) 

"* iot 

<J) (x, y, z, t) = <j> (x, y, z) e (44b) 


eq. (40) becomes 

2 - - ” 2Mico 7 co 2 “ _ . 

(l-M 2 )$ xil -f4. + r— + — <> = 0 (45) 

oo a 

oo 

Equivalently, eq. (45) results from applying a Fourier transformation to the time variable in 
eq. (43). Equations (43) and (45) present the governing equations for unsteady linear 
aerodynamics in the time domain and the frequency domain. The time domain approach has been 
valuable in giving insight into the nature of solutions for small values of nondimensional time. 
On the other hand, the frequency domain approach has provided the majority of the working 
methods used in aeroelastic design and analysis. 

The time domain form of eq. (43) is the basis of the "acoustic planform" analogy used in ref. 
97 to calculate the initial transient pressures and loads on airfoils and wings undergoing 
sinking and pitching motions. Figure 1 1 shows pressures on an impulisvely started sinking and 
pitching airfoil for several times during the first four chordlengths of motion. The Mach 
number is 0.8. The discontinuity in slope of the pressure indicates the progress of the 
upstream traveling pressure pulse generated at the trailing edge by the impulsive motion. This 
pulse travels at a speed of (a~ - U). The faster downstream traveling pulse (with velocity a~ + 
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U) can be seen in fig. 11(a). The resulting lift and moment transients for several Mach 
numbers are shown for the first 10 chordlengths of travel between 0 and 2 in fig. 12. 
Reference 88 discusses the use of piston theory in determining the starting value of pressure 
and integrated loads. At t = 0 the pressure is a function only of the local normal velocity of the 
surface, uo 


n - n = n a u 


/ A (* \ 

\ •+ o ; 


This leads to equations such as the following for the initial lift coefficient per unit angle of 
attack for any piane wing at any Mach number 


dc,(o) ^ 4 
da M 


(47) 


Reference 97 also contains time-accurate solutions of eq. (43) for 3-D wings in supersonic 
flow. 

In the alternative frequency-domain approach, extensive use has been made of fundamental 
solutions (Green's functions) of the governing equation such as the source and vortex potential 
functions for incompressible flow and the acoustic source pulse and doublet for compressible 
flow. Due to the linearity of eq. (43), superimposed distributions of these fundamental 
solutions are made to satisfy the boundary conditions. The assumption of simple harmonic 
motion, eqs. (44), enables the manipulation of the resulting expressions into functions which 
may be used to computed the strength of the singularity distribution. . 

For incompressible flow eq. (43) reduces to Laplace's equation 


V<j)=0 


(48) 


For two-dimensional flow and assumed harmonic airfoil motion, Theodorsen89 derived an 

analytical solution of (48). Garrick and Rubinow^O give a closed form solution of (43) for this 
same problem for supersonic flow while Possio's integral equation (see ref. 88) results from 
operating upon eq. (43) to produce a singular integral equation relating known downwash 
velocities, wa(x), on the airfoil to the unknown pressure difference, Apa (eq. 6-111, ref. 88) 
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w a (x) = - 


J Ap a (^) K (M, 

P JJ , b 


k(x-5) 

b 




- b < x < b 


(49) 


The kernel function, K, is composed of Hankel or Bessel functions and is a function of the Mach 
number and the assumed reduced frequency of oscillation, k. 

The corresponding singular integral equation for three-dimensional flow has been studied 
extensively (see, for instance, refs. 88, 91, 92) 

w(x, y) = — f f Ap a (^, q) K [M, k, (x - ^), (y - i))]d(;dTi ( 5 0 ) 

4n J S J 

Inversion of eq. (50) via substitution of assumed series expansions for the unknown pressure 
and numerical quadrature or collocation solution procedures have been used extensively for 
aeroelasticians. Knowledge of the functional behavior of the pressure loading near surface edges 
and slope discontinuities has been of great help in constructing suitable loading functions. 

Tables 4 and 5, from Ashley’s 93 survey, illustrate the singular behavior of ACp near wing edges 
and control surface boundaries. Note particularly the singular behavior for subsonic flow for 
the cases of control surface edges, (logarithmic, ~ In |x - xc|) and for wing leading edges (~ (x 

- xi)-i/2) and side/trailing edges (~ (y - yt) 1/2 . (xt - x)l/2). Aeroelastic forcing functions 
are calculated as weighted surface integrals of these pressure loading functions. Attempts to 
perform these integrations numerically without acknowledging the possibility of such singular 
behavior can lead to significant errors and influence solution convergence. Rowe and his 
coauthors 94 have approached this issue by modeling the singular behavior explicitly in their 
RHOIV "kernel function" computer code. With the singular portion evaluated analytically, the 
remaining integral equation is regular and straightforward numerical solution possible. An 

alternative method which is widely used is the doublet-lattice method of Albano and Rodden 95 . 
Here the lifting surface is divided into small trapezoidal panels (fig. 13). Within each panel 
line segments of acceleration potential doublets are placed on the panel quarter-chord line. The 
unknown doublet line strength for each panel is determined by satisfying the known downwash 
velocity boundary condition at the mid-point on each panel's three-quarter chord line. Thus the 
problem is reduced to a linear set of algebraic equations for the doublet line strengths. The 
choice of the 3/4 chord location for downwash evaluation is selected empirically upon noting 
that this results in the Kutta condition being satisfied. This paneling method can be easily 
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extended to multiple lifting surfaces as well as to body interference problems. 


Comparisons 96 of a convergence study comparing the RHOIV kernel function method and the 
doublet-lattice method are shown in figs. 14 and 15. Figure 14 shows the real and imaginary 
parts of unsteady lift and moment for an aspect ratio 4 wing pitching about its trailing edge at k 
= 0.3, M = 0.9. For this case 6 terms were required for the assumed pressure loading function 
in RHOIV in order to obtain convergence. The doublet-lattice results indicate that as usual the 
lift converges faster than the moment and that, for this case, the generalized force magnitudes 
converge quickly while the phase angles require more terms for convergence. Figure 15 shows 
a similar comparison for hinge moment due to control surface oscillations of a rectangular wing 


OA nffA I r» I I rf <*1 A A AAnfi/11 I K A ♦' A A TU /A AA oaA t. A A iJ A A | 4 A .... i . _ 

wiiiiwi ouiiauc uwimyuiauuii. iiio tswnuiuuio cut* iv = u.o, ivi = aiiu iu pressure terms were 

required for convergence in RHOIV. For up to 18 panels per chord ( and up to 6 panels on the 
control surface) the doublet-lattice result is converging slowly to the RHOIV result. For this 
case, the kernel function is more efficient in terms of computer resource units (CRUs). 

The behavior of solutions of (43) for Mach numbers near unity and/or for high reduced 
frequencies is of interest for the insight which may be gained of the transition between subsonic 
and transonic flows. Figure 16, from ref. 98, shows the pressure distribution on an airfoil 
oscillating in plunge for three values of k. The oscillations, which are most apparent for k = 5, 
are termed Kutta waves since, for this subsonic condition, they are generated at the trailing 
edge. For 2-D flow, isolated source solutions of eq. (43) can be viewed as cylindrical pressure 
waves radiating outwards from the source point at the acoustic speed, a~. Viewed from the 


translating airfoil, two wave fronts traveling at relative speeds of a~ ± U are seen on the z = 0 
plane. For high subsonic speeds and low supersonic speeds, the upstream propagating wave, 
a ~ - U, generates the surface pressure waves seen in fig. 16. The oscillations are in 

quadrature, with the real part having a cos [R (x - 1)] dependence and the imaginary part 
having a sin [R (x - 1)] dependence indicating an upstream propagating signal with wave 
number 


R = Mk/(1 - M) (51) 

radians per semichord. For M = 0.7 and k = 5 the wave number is 11.6 radians per semichord, 
agreeing with the oscillations shown in fig. 16. 

These pressure oscillations are also observed in solutions of (43) for low supersonic Mach 
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numbers. The pressure difference for oscillatory plunging 98 is 


r 8k 2 x 

AC d (x, k) = - [ 

p p 


_ ,, , — . . 4k -iwx 

f (M, xfo) + x — e 

° P 


T f^Xl h ° - 


lCOt 


152) 


Figure 17, from ref. 98, gives results for M = 1.02 and k = 0.4. In this case the oscillations 
are due to the overtaking of upstream traveling waves generated near the leading edge, leading to 
the supersonic wave number R = Mk/(M - 1) of 20.4 radians per semichord seen in the 
imaginary part over the forward part of the airfoil. Note the transition near midchord to an 
oscillation of one-half this wavelength caused by the interaction of the two terms of eq. (52). 
These pressure oscillations influence the integrated airfoil loads for Mach numbers near 

unity as shown in fig. 18 98 . The influence appears to be largest for low supersonic speeds. Also 
noted the reversal of trend of Re(ct) for M - 0.7 - 0.8, a feature important for aeroelastic 
analysis. 

In concluding this section, the equivalence of the time-domain approach of eq. (40) and the 

frequency-domain approach of eq. (45) is noted. Edwards 98 - 99 discusses several 
misconceptions regarding solutions of eq. (45) for diverging and converging oscillatory 
motions. Valid solutions of eq. (4) for arbitrary motions are obtained via Laplace • 
transformation of eq. (43). Solutions are functions of the Mach number and the generalized 
reduced frequency § = h/U (a + ico). Generalized unsteady aerodynamic methods are becoming 
widely used for aeroelastic analysis and in the design of active aeroelastic control systems. 

Figure 19 98 shows a typical result, with the unsteady lift coefficient of a wing given for rigid 
plunging motions at M = 0.5. Results are given as a function of amplitude, r, of motion for 

three phase angles, 9, where s = re' 0 . Harmonic oscillation results are given by 9 = 90° while 
converging and diverging motion results are given by 120° and 60° respectively. 


Accuracy and Resource Requirements 

The previous section briefly surveyed solution methods which are available for the linear 
potential equation. These methods are well developed and, for the most part, do not tax available 
computational resources. This is definitely not the case for solution methods of the nonlinear 
fluid dynamic models to be considered in the remainder of these lectures. There are, in general, 
no closed form solutions available for these equation sets and iterative numerical methods are 
the rule. Computational unsteady aerodynamics requires careful selection of equation level 
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based upon considerations of the flow physics involved, required accuracy and the number of 
cases required to perform aeroelastic analysis. The discussion of figures 1-8 has summarized 
the range of flow physics involved in current aeroelastic problem areas. Accuracy 
requirements are dependent upon the type of flow as delineated in fig. 3. Table 6 summarizes 
the current accuracy of predictions of key aeroelastic response mechanisms and suggests 
accuracy requirements for computational aeroeiasiic analysis. Note that the current accuracy 
level for high subsonic speed attached flow, the most relevant gage of current capability, is on 

the order of 10 percent. Novel computational methods will have to do significantly better than 
this to be competitive. 

The resources required are a function of the number of computer operations required per 
case and the number of cases to be calculated. Tables 5 and 6. from ref. 100, summarize the 
resources required for typical aeroelastic analyses. 

Table 7 indicates the computer resources required to perform a flutter analysis of a 
complete aircraft configuration at one Mach number. Time-marching transient aeroelastic 
response calculations are used to determine the flutter condition. This involves, on average, 
four response calculations: two to calculate steady flow field conditions and two transient 
responses bracketing the flutter speed. Modal frequency and damping estimates from the 
responses are determined and the flutter speed interpolated from the damping estimates. 
Calculations have been performed for a complete aircraft configuration with a transonic small 
disturbance (TSD) potential code using 750,000 grid points. The calculation of one flutter 
point for this case on the CDC VPS-32 computer would require 2.3 CPU hours. Estimates of 
similar calculations using the full Navier-Stokes equations would require 77.8 CPU hours. 
Conditions for this estimate are; a Reynolds number of 10 million, 7 million grid points and an 
assumed computational speed of 100 million floating point operations per second (MFLOPS), 
Table 8 summarizes computational requirements for flutter calculation of a complete flutter 
boundary for wing/body/canard configuration on the CDC VPS-32 computer operating at 100 
MFLOPS and on the NAS CRAY II computer operating at 250 MFLOPS. Again, four response 
calculations per flutter point are assumed. It is assumed that ten flutter points will be 
calculated to define the flutter boundary versus Mach number. The left hand column indicates 
the difficulty of the flowfield calculation as defined in figure 1; type I for attached flows, type II 
for mixed (alternately separated and attached) flows and type III for fully separated flows. The 
second column indicates the fluid dynamic equation level needed to accurately model the flow 
physics of the problem. Note that two-dimensional strip boundary layer models are assumed 
for interactive viscous-inviscid calculations for the potential and Euler equation methods. It is 
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anticipated that potential equation models will be adequate for flutter calculations of type I 
attached flow conditions and may also be quite useful for some type II mixed flow cases. Full 
potential equation codes will require about 50 percent more computer resources than TSD 
methods due to the necessity of conforming, moving grids, among other considerations. Euler 
equation methods should also be adequate for these conditions and, in addition, be able to treat 
more difficult type III fully separated flows. Euler equation methods are estimated to require 
approximately twice the resources of TSD methods. The full Navier-Stokes equations, which 
should only be required for type II and III flows, require approximately 30 times the resources 
of the Euler equations (at a Reynolds number of 100 million). 


29 


PART III COMPUTATIONAL SOLUTION METHODS 

In this lecture, methods which have been developed for the solution of the fluid dynamic flow 
models described in Part II will be discussed. The presentation will be partly historical, 
leading from the simpler TSD potential equation to the higher equation level models. Although 
not always the case, this sequence broadly reflects the progression of research in this area. 
Solution algorithms will be discussed and results both of historical interest and indicative of the 
current state will be presented. It is helpful to regard the evolution of progress in this area as 
following four broad stages: 

i . ) Early computational demonstrations 

ii. ) Maturation of computationai methods 

iii. ) Application to realistic configurations 

iv. ) Type II mixed flow computation 

The first three stages are those typically encountered in the development of any novel technology 
within a given problem area. Jn the current context, they relate to development of capability of 
type I attached flow (fig. 3 and Table 4) which has been a dominant focal point of applications. 
The fourth stage listed serves as a caution to too rigid a categroization as parallel efforts for 
type II mixed flow calculations have proceeded along side, and sometimes ahead of those for 
attached flows. 

A thread which may be discerned in reviewing this field is the continued evolution of 
computational methods, with applications and evaluations by comparison with experiment, to 
successively more difficult cases. Thus, for instance, an algorithm introduced to treat type I 
attached flow cases is upgraded in capability to enable treatment of more difficult type I and 
type II cases and possibly even some type III cases. 

In the following, results of calculations for a number of the "Computational Test” (CT) cases 

drawn from the AGARD Standard Aeroelastic 2 -D 16 and 3-Di? configurations will be presented. 
These experimental cases are mostly for harmonic oscillations of a wind tunnel model in a rigid 
(i.e. the models have been made as stiff as is practical in order to minimize aeroelastic 
deformations) degree-of-freedom about a steady mean condition which, for pitching 
oscillations, is described as 

a(t) = a m + a o sin kt (53) 

The reduced frequency is based upon reference semichord unless otherwise specified. 
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Response quantities, such as lift coefficient, are measured or calculated for a sufficient 
number of cycles such that a stable average for the time history of the response is available. 
The response will, in most cases be cyclical at the frequency of the forced motion, k, but it will 
in general not be purely harmonic. That is Fourier analysis of the averaged response of, for 
example, the lift coefficient will give 

Cj(t) = Cj + ^ [Re n (Cj) sin nkt + Im n (Cj) cos nkt]; 0 < kt < 2 jc (54) 

° n=l 

with both Re n (Cl) and lm n (Cl) having nonzero values for n > 1 . Nevertheless, the fundamental 
response at the forcing frequency, given by Rei (Cl) and Imi (Cl) is most important since it is 
most easily measured and almost always accounts for the largest part of the response. The 
higher harmonics in the response given by the Fourier coefficients with n > 1 indicate the 
degree of nonlinearity in the aerodynamic response. 

The Fourier coefficients in eq. (54) may be obtained by several different approaches. If a 
number of the harmonic coefficients are desired in order to study nonlinear behavior, then Fast 
Fourier Transform (FFT) methods are most efficient. On the other hand, if only a few 
coefficients are of interest, the traditional evaluation via numerical integration is efficient. 
Finally, the following simple estimate of the fundamental response is useful 

Re i ( C i> * T [ C i ( l ) - C, (t) 3 ]-C (55a) 

L kt = -£- kt = — n 

2 2 


Im 1 (C 1 )« I [C 1 (o)-C,(t) - C, 


(55b) 


TSD Potential Equation, 2-D 

The two-dimensional TSD potential equation has been extensively studied at several levels of 
approximation. Reference 101 gives details of an alternating direction implicit (ADI) solution 
algorithm for the equation 


c <J) tl + A <> xt = B <>xx + 


ZZ 


(56) 
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The potential is normalized by cUS 2/3 where c is the airfoil chord, and 8 is the airfoil thickness 
ratio. The coordinates x, z and time, t, are normalized by c, c/8 1/3 and or 1 where w is the 
frequency of unsteady motion. The coefficient A = 2 k M“78 where kc = coc/U. In their 

LTRAN2 code, Ballhaus and Goorjian^ implemented a solution of the low frequency version of 
(50) by setting C = 0 and eliminating the time derivative terms the boundary conditions, eqs. 
(37,38). In ref. 14 B = (l-M)/5 - M m (Y + 1) <t> where the choice of the exponent is 

oo oo ^ 

arbitrary. Reference 14 made m a function of M- such that the critical pressure 

* * 
coefficient, C p , predicted by eq. (56) matched the exact isentropic C p . Houwink and van der 

Vooren 102 extended the range of applicability of LTRAN2 by adding the time dependent boundary 
condition terms, and solving a modified TSD equation wherein B = (1 - M )/5 - M (y + 1) 6 

OO OO A 

where y* - 2 - (2 -y)M 2 . The resulting code was termed LTRAN2-NLR. Whitlow 101 

oo 

extended the LTRAN2-NLR code by implementing Rizzetta and Chin's 103 solution of the complete 
TSD eq. (50) where C = k 2 M 2 /5 2/3 . The resulting code was termed XTRAN2L 104 . 

C oo 

The ADI method advances solutions from time step n to time step n + 1 using successive 
sweeps in the x and z directions 


x-sweep: 


— 5 (6. . - 6 n .) = D f. . + 5 <t> n . 

^ x VY t,j Y i.j / x ij zz Y i,j 


(57a) 


z-sweep: 



where <|> is an intermediate level potential and 


56 ='■ 


X. , 
l+l 




NY i,j Y i-1, y 
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(57b) 
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, 1 - NT 

f . . = - [ — — ^ <b n . + B n . A ] 
•o 2 1 S 2/3 Y,,J ,,J Vx iJ 


1-M „ . 


The mixed difference operator, Dx, is constructed to maintain conservation form. Murman-Cole 
spatial differencing results in the following form for Dxfi.j 


D f .. = 


x x. , - x. , 

i+i i-i 


[(1 -A)<WW +e i-i ( f i-i/2.j - w ] 


(58) 


£ i = 


0 


tf+m.j + Qm.j > ® 


^+l/2,j + ^i-l/2.j < 0 


With -the LTRAN2 code, Ballhaus and Goorjian 14 were able to reproduce Tijdeman's 3 type A, 
B and C shock motions for the NACA 64A006 airfoil with an oscillating flap. The motions were 

also demonstrated computationally by Magnus and Yoshihara 13 and Magnus 7 using an explicit 
Euler equation code. Figure 20, from ref. 14, gives these three calculations. The computational 
conditions for these cases are 


Iy os. 

Mach 

K 

flap amplitude 

A 

0.875 

0.234 

1.0 deg. 

B 

0.854 

0.179 

1 .0 deg. 

C 

0.822 

0.248 

1.5 deg. 


These conditions are 0.15 - 0.28 Mach lower than Tijdeman's test conditions, very likely due to 
wind tunnel wall interference. The Euler code used a Lax-Wendroff explicit differencing 
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solution algorithm and a Cartesian grid with an embedded fine mesh around the moving shocks 
was used. The boundary conditions were applied at the mean airfoil position. The Euler code 
used 5484 grid points and required 1500 seconds per cycle whereas the LTRAN2 code required 

8 seconds (CDC 7600 computer)i4. This significant reduction brought the expense of 2-D 
unsteady transonic CFD calculations within the reach of many researchers. 




V 


Monotone Differencing . It has 


shown that Iv'urman-Coie differencing allows stable 


entropy-violating expansion shocks to be computed as part of the numerical solution^. 

Reference 104 showed that it can also trigger numerical instabilities. Figure 21101 shows 
such a case for the MBB A-3 airfoil oscillating in pitch about its leading edge at M = 0.8 and 
k = 0.2. At kAt = 254° an instability is developing on the lower surface at the leading edge that 
causes program failure within several iterations. When the monotone differencing scheme of 

Engquist and Osherioe j s used, expansion shocks are not admitted and significant increases in 
allowable time steps over those allowed by the Murman-Cole scheme are achieved. The Engquisl 
Osher (E-O) scheme was first used in implicit algorithms by Goorjian et al.105 anc j a S j m j| ar 
implementation is described by Whitlowioi. The E-0 method is incorporated into the ADI 
procedure by modifying the mixed difference operator D x in eq. (58) to: 


D x f i- l/2.j A x h-m.i + A x f i- 1/2. j 


i+1 


i-l 


(f . f 

Vl i+l/2.j i-l/2,j ' *i-l/2,j 


A 

+ f. 


' f i-3/2.p 


(5 9) 


The f and f operators are given in ref. 1 01 . 

In the following, examples of results from the XTRAN2L code will be discussed and 
additional modifications to the ADI algorithm will be described. The nonlinear term, <t> x <D XX , 
leads to the formation of shock waves for transonic speeds which the numerical differencing 
solution schemes are designed to capture. Pressures in the vicinity of a shock can vary in a 
quite nonlinear fashion with respect to surface motion and the importance of this nonlinearity 
in aeroelastic applications is of great interest. Typically, applications of unsteady 
aerodynamics involve not the local unsteady pressures but their integrated value, weighted by 
shape functions describing surface motion. In 2-D the lift and moment coefficients resulting 
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from airfoil pitching and plunging motions are the relevant airloads. It has been observed that 
these airloads behave in a sufficiently "locally linear" manner that linear harmonic analysis 
methods are quite useful. 

Pulse Transient Unsteady Airload Calculations . If the airload response to surface 
motion is sensibly locally linear then airload frequency response functions may be determined 
from Fourier transform analysis of transient responses. References 104 and 107 describe the 
pulse transient method which is based upon this assumption. Starting from a converged steady 
state solution the airfoil boundary condition eq. (37) is prescribed to simulate an exponentially 
shaped pulse motion (e.g. pitch). 

a(t) = a o + ctj e ' w(t ' ^ (60) 

Fast Fourier transforms (FFTs) of the lift and moment coefficients and the angle of attack time 
histories are calculated. The lift and moment FFTs are divided by the angle of attack FFT to 

obtain the do and cma frequency response functions. Figure 22io 0 shows a typical pulse 

transient result for the NACA 0010 airfoil at M = 0.78 and ao = 0.0°. The figure shows only 
the early portion of the transient to illustrate the fluid dynamics resulting from the pulse. 
Calculations are typically continued for ~ 1000 time steps to allow the loads to return to steady 
state. Figure 23 shows the resulting airload transfer functions for this case and for similar 
thickness NACA 64A010 and parabolic arc airfoils. This figure indicates the degree of 
difference in unsteady airloads which results from the nonlinear transonic steady flow 
condition. Figure 24 indicates the correlation between airloads obtained using this method and 
the harmonic oscillation method for a six percent parabolic arc airfcil at M = 0.85. The pulse 
shape is chosen to give reasonable results for reduced frequencies up to k ~ 2.0. Its use gives 
considerable detail in the frequency domain from a single transient calculation resulting in a 
considerable reduction in cost over the harmonic oscillation method. 

In order to obtain accurate results with the pulse method careful attention to numerical 
details is required. In order to avoid spurious low frequency results, the starting value of a(t) 
in eq. (54) should be equal to ao. For realistic pulse amplitudes, ai, w and t c should be chosen 

so that a (o) - ao = ai e- wt c**2 < 10- 6 . Also, the method is quite sensitive to lack of 
convergence in the starting steady-state solution. Drift in the unforced (ai = 0) lift coefficient 
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of 10-3 - 10-4 C an again cause significant low frequency errors. 


Computational Grid Dyn amic Effects 

The pulse trnasient method has been an effective tool in studying the effect of the 
computational grid upon unsteady calculations. Reference 107 uses this pulse transform 
technique to demonstrate key features of the relation between computational grids, boundary 


conditions, and dynamic calculations. The importance of controlling reflections of disturbances 
from the outer grid boundaries and from internal grid points is shown in fig. 25. In order to 
compare with linear theory, the case of a flat plate airfoil at M = 0.85 is shown. Three lift 
coefficient responses resulting from quickly pitching the airfoil from 0 to 1 degree and back to 


0 are given. In figs. 25(a) and 25(c) the default XTRAN2L grid given above was used while in 
figure 25(b) an exponentially stretched grid extending ±200c in x and ±2327c in z was used. 
The latter grid contained 113 x 97 points in the x and z directions. The default XTRAN2L 


gridi 07 is 80 x 61 points in x, z and covers a fixed physical domain extent of ±20c in x and 
±25c in z. An algebraic grid stretching is used to distribute grid points between the airfoil and 
the outer boundaries. This grid point distribution was selected to alleviate disturbances which 
can be generated in regions of large grid stretching. On the airfoil the x-grid has 51 grid points 
having a uniform spacing of 0.02c with an additional point near the leading edge. Figure 25(c) 
was obtained using the non-reflecting boundary conditions given by Table 3 while figs. 25(a) 
and 25(b) utilized reflecting boundary conditions. Of particular importance are the outer 
z-boundaries. The disturbance at x - 40 in fig. 25(a) correlates with the acoustic propagation 
time for travel to and return from these boundaries. The option of moving these boundaries to 
large distances, as in fig. 25(b), introduces the complication of severe grid-stretching in the 
near-field. In this case, reflections from the outer boundaries do not occur, but disturbances 
seen from x « 20 to 50 correlate with propagation times for travel to and return from regions 
of the z-grid where grid spacing first becomes more than two chordlengths. Neither of these 


anomalies is seen in fig. 25(c). 

Figure 26 gives the clot frequency responses calculated from these transient responses. 


Reflections from the outer z-boundary, fig. 26(a), contaminate the unsteady airloads at low 
reduced frequencies, k < 0.15, whereas the distrubances originating from the near-field grid 
stretching, fig. 26(b), contaminate the airloads in the frequency range 0.2 < k < 1.0. Figure 
26(c) shows that excellent agreement with linear theory can be achieved for moderate 
frequencies. Other calculations verify that these features, which are most easily studied for 



36 


linearized examples, carry over to nonlinear transonic calculations. Reference 109 gives 
further examples of computational grid effects for unsteady 3-D calculations. 


Viscous-lnviscid Interaction . The inviscid TSD equation (56) does not incorporate 
viscous effects which can be important for high speed flows. It is possible to account for 
unsteady viscous effects by coupling a viscous boundary-layer model with an otherwise inviscid 
analysis. As commonly implemented, the inviscid outer flow solution provides the surface 
pressure distribution needed to solve the boundary layer equations. This yields the boundary- 
layer displacement thickness distribution which is used to modify the airfoil surface tangency 
boundary condition for the next iteration of the outer inviscid flow solution. 

Howlett 10 describes such a method implemented in the XTRAN2L code. The effect of a viscous 
boundary layer for attached turbulent flow is modeled in a quasi-steady manner by means of 

Green’s lag-entrainment equations 1 1 1 as implemented by Rizzetta 112 . In this integral method 
the displacement thickness 5* is computed as a function of the boundary-layer momentum 
thickness 0 and the shape factor H: 

6* = 0 • H (61) 

The functions 9 and H are determined together with the entrainment coefficient Ce from Green's 
lag-entrainment equations. In the nondimensional variables consistent with eq. (56) these 
equations are 


.£(!) = f + f<j> =_l-(H + 2-M 2 )e 2/3 ^4 

dx v c / 1 2 r xx 9 v e' r Y xx 


(62) 


f sr 1 - <3 + f A* = < c e • -r c t) is: + H . (H + » f ♦» 


( 63 ) 
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+ <ir^rW- F[1+a075M J 

e 


1 +Li r M 2 

T ]e ™^ 

1+0.1M c 


XX 


(64) 


The subscript e in these equations refers to the quantities evaluated at the boundary-layer edge, 
the subscript EQ denotes the equilibrium conditions, and the subscript EQO denotes the 
equilibrium conditions in the absence of secondary influences on the turbulence structure. 
Expressions for the functions contained in eqs. (62) - (64) are listed in Appendix A of ref. 
100 . 


n/MimrtrAnm r\f fhn troilinn.orlno tho como oniiotinnc 3rci annlioH t n oarh ciHo nf tho \a/pUp 
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surface independently with the skin friction set to zero and the dissipation scale length doubled 
to account for the observed far field behavior of wakes. 

Coupling between the boundary-layer and inviscid analysis is through the boundary 
conditions on the airfoil and wake, eqs. (37) and (39), which are modified to 




(65) 


Mg = [8*/5c) x ] (66) 

Coupling between the inviscid analysis and the boundary-layer is through the quasi-steady 
pressure gradient, <J>xx, in eqs. (62) - (64). Explicit coupling between the boundary layer and 
the inviscid solution is used for the airfoil boundary condition, eq. (65), since this allows a 
substantial increase in the allowable time step. That is, the last term in eq. (65) is evaluated 
as 


i O. - - O. - 

(S' /6 c) .X-SL-bL 


8c 


i+i 


x. 


i-l 


(67) 


In the next section results calculated with the XTRAN2L TSD code are compared with 
experimental results. Results from the interactive viscous-inviscid model will be designated as 
IV-TSD 
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Two-Dimensional Transonic Calculations 
In this section, comparisons of calculated and experimental unsteady pressures are given for 
the NACA 64A006 airfoil, the NACA 64A010A airfoil and the NACA 0012 airfoil. The cases are 

chosen from the AGARD Standard Configuration Computational Test (CT) cases 18 presented in 
refs. 100 and 113. 

NACA 64A006 . Tests of this airfoil 8 all involve oscillation about zero mean angle of a flap 
with hinge axis located at three-quarter chord. The steady flow pressure distributions for Mach 
numbers from 0.80 to 0.875 are shown in fig. 27. At M = 0.85 the IV-TSD model results agree 
well with the data while at M = 0.875 the viscous results correct roughly one-half of the 
discrepancy in the inviscid shock location. Post-shock pressure levels are well predicted by 
the IV-TSD model. Figure 28 shows the displacement thickness, 5* for these cases. The 
thickening of the boundary layer by the shock is apparent for M = 0.85 and 0.875. Figure 29 
shows unsteady upper surface pressures at a low reduced frequency, k = 0.06 and at a moderate 
reduced frequency, k = 0.24. Figures 30 and 31 show the corresponding airloads for these 
cases. In general, the agreement between experiment and calculations improves with decreasing 
Mach numbers and increasing frequency. As for the steady pressures, the unsteady results in 
fig. 29 indicate that the IV-TSD model accounts for a portion of the discrepancy in the shock 
pulse. The airloads show reasonable agreement with the data for the lower Mach numbers with 
trends also reasonably predicted. However, with the onset of significant transonic effects at 
M = 0.85 it is obvious that further computational improvements are called for. Note, in 
particular, that linear theory continues to perform well for these cases. 

NACA 64A010A . These cases are for the model tested at the NASA Ames Research Center 62 in 
which the model had a small amount of camber and was 10.6 percent thick. These cases are for 
the model pitching about the quarter-chord at nominally zero degrees pitch angle. Figure 32 
gives the calculated and measured steady pressures for M = 0.5 and 0.796 and the displacement 
thickness for these cases is given in fig. 33. At the lower Mach number, agreement is very good 
with almost no viscous effect evident. At the higher Mach number, the viscous results move the 
shock forward several percent chord and slightly weaken the shock strength. The data sets 
encompass Reynolds numbers variations ranging from 2.5 x 10 6 to 12.5 x 10 6 based upon 
chord and fig. 33 shows the resulting predicted variation in displacement thickness. Figure 34 
shows the upper surface unsteady pressure for M = 0.5 and k = 0.10. There is excellent 
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agreement for this case which is typical of other results for this Mach number. Figure 35 
shows upper and lower surface unsteady pressures at M = 0.796 for reduced frequencies from 
0.025 to 0.3. Pressure levels ahead of the shock are generally well predicted. The inviscid 
shock pulse is too strong and too far aft, with both effects being improved by the IV-TSD model. 
The viscous model also improves the post-shock pressure comparisons. Still, there is again an 
evident need for further increases in accuracy in the region of the shock which is also apparent 
in the unsteady airloads shown in fig. 36. Note the characteristic overprediction of (cma) at 
low values of k only a portion of which is corrected by the viscous model. Similarly, the values 
of Re (cma) at low values of k appear to have an anomalous trend. Trends due to amplitude of 

oscillation are well accounted for as is shown in fig. 37. For amplitudes of 0.5, 1 .0 and 2.0 
degrees, the IV-TSD model very nicely corrects the inviscid shock pulse signature as well as the 
post-shock pressures. 

NACA 0012 . This case involves a 1 2 percent thick symmetrical airfoil tested with free 
transition for sizable mean angles and oscillation amplitudes as well as cases with constant pitch 
rate ramping motions to high angles, ref. 18. Results for the latter cases are shown in ref. 113 
and demonstrate that the TSD code yields surprisingly good lift coefficient estimates for these 
transient ramping motions up to angles near stall, a - 8-10°. Figure 38 presents results for 
the former cases, with total lift and pitching moment coefficients plotted versus pitch angle. 

The first three cases are for oscillations of 2.5 and 5 degrees about non-zero mean angles while 
the last case is for oscillations of 2.5 degrees about a zero mean angle. Agreement for the lift 
coefficients varies from good to very good whereas the moment coefficients for the first three 
cases show a systematic difference with experiments due to underprediction of pressures near 
the leading edge suction peak in pressure. The characteristic shape of the cm - a curves for the 
first three cases is due to a large second harmonic contribution. In contrast, the shape in the 
fourth case is due to an increased third harmonic component. Viscous results from the IV-TSD 
model are shown for the first and fourth cases. In the first case, viscous effects produce what 
appears to be a decreased agreement with experiment. The second and third cases were not 
amenable to the "direct" IV-TSD model described above since they involve incipient separation. 
Such cases require an "inverse” boundary layer method to handle mildly separated flows. 

Parametric Studies . A number of parametric studies of transonic 2-D unsteady 
aerodynamics have been published. Bland and Edwards 114 investigated the effects of airfoil 
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shape for the 10.6 percent thick NACA 64A006 airfoil and the 8.9 percent thick MBB A-3 
airfoil. The latter is a supercritical airfoil with significant aft loading. Figure 39 shows the 
steady pressure distributions on these airfoils for a range of Mach numbers and fig. 40 shows 
the unsteady pressure magnitude and phase angle due to pitching for the same Mach numbers fcr 
k = 0.15. It was observed that the frequency response functions of these two airfoils were very 
smilar when a Mach number shift of 0.01 was used to account for the thickness difference. 

Figure 41 shows this comparison for airfoil pitching. Using this Mach number shift, it was 
shown that the two airfoils exhibited similar transonic flutter characteristics for structural 
dynamic parameters representative of a swept-back wing section. 

The effect of amplitude upon transonic unsteady aerodynamics was also studied 1 14 . Figure 
42 shows pressure distributions due to pitching for the NACA 64A010A airfoil at M = 0.7 and 
k = 0.15. For oscillation amplitudes from 0.25 to 2.0° normalized pressures ahead of and 
behind the shock are little affected by varying amplitude while the shock pulse is smeared out 
with increasing amplitude. Note, however, that the area under the shock pulse remains nearly 

constant, lending further credibility to the use of locally linear analysis methods. Howlett 110 
also investigated amplitude effects with his interactive viscous boundary-layer method and 
found no significant nonlinear effect for a moderate transonic case. Figure 43 shows lift due to 
pitch for the NACA 64A010A airfoil at M = 0.796 for 0.1 and 4.0 deg. pitching amplitude. 

Batina 108 investigated the effects of airfoil shape, thickness and angle-of-attack using the 
XTRAN2L code. As a reference, fig. 44 gives the linear cma frequency response function at M = 

0.80 (linearized aerodynamics and flat plate airfoil). Shape effects were investigated using 
three 10 percent thick airfoils; the NACA 0010, the NACA 64A010 and a parabolic arc airfoil. 
Figure 45 shows steady pressures for M = 0.76, 0.78 and 0.80. Note the formation of shocks 
near the locations of maximum thicknesses at 30, 40 and 50 percent chord. Figure 46 gives the 
Cma frequency responses for these cases. Comparing figs. 44 and 46 two transonic features are 

apparent. First, the development of nonzero values for Re (cma) at k = 0 with increasing Mach 

number reflects the aft shift of the steady center of pressure for the NACA 64 A0 10 and 
parabolic arc airfoils. Secondly, the "wave" feature, seen most prominently in both the real and 
imaginary responses near k = 0.6 for M = 0.78, is novel. It is apparently independent of 
airfoil shape and occurs at lower reduced frequency for higher Mach number. 

Reference 108 also studied airfoil thickness effects. In steady flow, we have a transonic 

similarly rule 83 whereby similar airfoils with equal similarity parameters 
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Figure 47 shows such scaled steady pressures at three similarity scaled Mach numbers for the 
NACA 0008, NACA 0010 and NACA 0012 airfoils. The similarity parameter values xi = 

1.4749, 1.3270 and 1.1851 lead to, for example, Mi = 0.76, 0.78 and 0.80 for the NACA 
0010 airfoil. The corresponding cma frequency responses are given in fig. 48. Although no 

similarity rule is available for the unsteady case, it is interesting to note that, for low and high 
frequencies, these scaled airfoils behave quite similarly. In contrast, the wave feature which 
will be identified with "aerodynamic resonance" below does not scale according to eq. (68). 
Reference 1 points out that such resonances are related to the shock motion types identified by 
Tijdeman. Figure 49 shows CmS frequency responses for the NACA 64A006 airfoil for the two 

Mach numbers at which type B and C shock motions have been calculated (fig. 19) for the 
frequencies denoted by the symbols. The type B motion corresponds to the resonance frequency 
of the aerodynamic resonance feature (the maximum amplitude of Im (Cm5) for this case) at 

M = 0.854 whereas the type C shock motion frequency is below the M = 0.822 resonance 
frequency. 

Finally, Edwards et al. 115 present parameter studies of 2-D transonic flutter 
characteristics. The effects of airfoil shape and angle-of-attack are noted and detailed results 
for a variety of structural dynamic integration methods are shown. Effects of motion amplitude 
and time step size are also noted. 


Shock Generated Entropy 


Thus far, transonic unsteady aerodynamic calculations have been presented which 
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demonstrate the capability of the TSD potential equation to capture nonlinear shock waves and to 
predict trends for transonic unsteady airloads. An evaluation of figures 29 - 38 indicates that 
such computations are not accurate enough to displace linear theory. Augmentations to the basic 
TSD equation and/or resorting to higher equation levels are called for. In this section a 
significant augmentation to the TSD equation is described. 

The conservative full potential and transonic small disturbance potential equations are 
derived with the assumptions that the flow is irrotational and isentropic. While it is recognized 
that entropy is generated within shock-waves, the use of potential theories to study transonic 
flows with weak shocks has progressed assuming that this entropy generation was a higher 
order effect. It is now understood that disregarding this effect can lead to serious disagreement 
with more exact solutions for physically interesting situations. 

A common approximation in formulating the full potential equations is to impose 
conservation of mass and energy while satisfying the isentropic and irrotationality 

conditions 116 . Reference 116 shows that the shock jump conditions for such a "conservative 
potential" equation deviate from the Rankine-Hugoniot shock conditions as the Mach number 
ahead of the shock increases. In ref. 116, the implications of this effect were studied by 

calculating ci versus a for a range of Mach numbers. It was known 117 that the symmetric NACA 
0012 airfoil at a » 0° exhibited multiple solutions for 0.82 < M < 0.85. Figure 50 shows that 
such ranges of multiple slutions can be found for all Mach numbers for sufficiently large angles 
of-attack. More importantly for transonic aeroelasticity, it is concluded that well before a 
reaches values at which multiple solutions occur, the lift-curve slope, cia, can become 
unphysically large. 

Williams et al. 118 have investigated the effect of nonunique solutions of the unsteady TSD 

equation. Figure 51 118 gives three different calculations of lift coefficient versus a for the 
NACA 0012 airfoil at M = 0.85. Figure 51(a) gives the upper surface pressure distributions 
for the three multiple solutions indicated by A, B, and C in fig. 51(b). Solution B is a 
symmetric nonlifting solution while the other two are lifting solutions. Figure 51(b) gives the 
lift coefficient versus angle-of-attack for 1.) Quasi-steady conditons, k = 0, and pitching 
oscillations for 2.) k = 0.01 and 3.) k = 0.05. Solution B is not a stable solution and diverges 
with an extremely large time constant to either A or C depending upon initial conditions. At k = 
0.05 the solution oscillates about the positive lifting solution. While the average lift curve 
slope is not unreasonable the solution must be regarded as anomalous. In contrast, the solution 
for k = 0.01 exhibits a hysteresis loop, jumping between the two stable steady solutions. The 
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large phase lage implied by this slution is unphysical and caution must be exercised against such 
calculations. 

Fuglsang and Williams' 1 19 implemented a nonisentropic formulation for the 2-D TSD 
equation. In the 2-D form of eq. (30) the streamwise flux may be written as 

2 . 1 _ 2 

t = fl - M )U - — HU ( / U ) 

Zr 

where u = <j>x. The coefficient B may be chosen in a variety of ways but should be asymptotic to 
(y + 1) as M approaches unity. Reference 119 replaced eq. (70) with 

2 * V 2 

f=(y+l)M 2 R(VV -— ) (71) 
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and asterisked quantities refer to sonic conditions. This new flux is identical to that of eq. (70) 

to 0(u 2 ) when u is small. The new form is derived from a formal asymptotic expansion of the 
Euler equations, including the effect of shock-generated entropy. A modified pressure 
coefficient to account for this entropy and a modified wake boundary condition to account for 
entropy convection complete the nonisentropic modeling. The pressure coefficient is modeled as 


( 72 ) 



where 
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V-2(V4\> 

is the linear isentropic term and Cps is the correction due to the entropy jump 


P ’ y(y+l)M 2 

The entropy jump is evaluated using the computed velocity upstream of the shock and the 
Rankine-Hugoniot normal shock jump condition 


(73) 
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where 


e = (Y+ 1)/(Y- 1) 

_ _ 2-2 
r = u 2 /Uj = R /Uj 

“ = 1+ V u , 

The wake boundary condition, eq. (38), is modified to account for the entropy as 

[4 ) x + <|) l ]=y[ c p j ] (76) 


The entropy is assumed to be convected down the wake at the freestream speed leading to 
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dt p« dx p* 


( 77 ) 


Figure 52 shows isentropic and nonisentropic results compared to Euler and full potential 
results for M = 0.84 and a = 0.25. The nonisentropic results are very similar to the Euler 
results. No mulitple solution conditions have been observed with the nonisentropic model and 
values of lift-curve slope are reasonable. Also, low frequency unsteady calculations do not 

exhibit the hysteresis effect shown in fig. 51(b) 119 . 

TSD Potential Equation. 3-D 

The success of the ADI solution algorithm in enabling efficient 2-D calculations led to its 

extension to 3-D in the XTRAN3S 28 - 29 co <j e . In this code, the physical grid lines in the x, y 
plane conform to the wing planform and the grid is extended in regions off of the wing planform 
using a shearing transformation to map onto a rectangular computational domain: 

% - $ (x, y). Tj = y, C = z> * = t (78) 

In computational space eq. (30) becomes 120 
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The time-accurate solution of eq. (79) via the ADI algorithm is summarized by Borland and 
Rizzetta 29 : 
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where 8^, 8t|, 8(£ are second-order accurate central difference operators, 8^ is a first-order 


backward difference operator, is a mixed-difference operator based on the 


sign of (a" + 2b 4>£) and Dti is a mixed-difference operator based on the sign of 
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Thus numerical stability is maintained using Murman-Cole differencing. Expressions for a, b, 
X and Y are given in ref. 29. This algorithm is used to advance the solution for <(> from time t to 


time t + At (i.e. 4» n to <|> n+1 ). All of the terms contributing to the streamwise portion of the 
equation are treated implicitly, as well as the second difference Sr> (c 8 r> 4 ») and c 8 t£ <t>- The 


remaining cross terms contained in the expressions for X and Y are handled in an explicit 
manner. Because of this there is an inherent time step limitation for stability of the 3-D 
method not present for fully implicit methods such as the 2-D ADI method of ref. 14. 

A comparison of calculations from the XTRAN3S code with experimental data from a 
rectangular supercritical wing oscillated in pitch is given in ref. 109. Figure 53(a) shows 
steady pressures for two span stations at M = 0.7 and a = 2°. For this low transonic condition, 
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the agreement with experiment is good except near the leading edge. The unsteady results in fig. 
53(b) are in good agreement with experiment over most of the chord. Inboard, there is an 
improvement in the prediction of the leading edge suction peak over the linear RHOIV result. 
Outboard, there is evidence of viscous effects in the phase angle in the aft cove region. 


K'vvmgq vjrig i rgngfQrmgtipq . The semi-implicit ADI algorithm used in the XTRAN3S code 
required increasingly small time steps to maintain stability as wing sweep and/or taper were 

increased. Guruswamy and Goorjian^i modified the original grid shearing transformation, eq. 
(20), which resulted in highly skewed outer grid boundaries for such cases. The modified 
transformation maps from a rectangular physical domain onto a rectangular computational 
domain while maintaining the alignment of the grid with the wing planform. Reference 120 
describes the details of a similar transformation method. In regions upstream and downstream 
of the leading- and trailing-edges and their extensions, the grid spacing is given by a cubic 
shearing function, which for the. upstream region is 


\ = X LE + A l i „ + A 3 V *u = **•* 


uP 


( 81 ) 


The coefficients Ai and A 3 are selectedi20 to enforce smooth grid metrics for the grid cells 
adjacent to the wing. Outboard of the wingtip the leading and trailing edges are extended using 
similar cubic functions that match the leading and trailing edge slopes at the tip and intersect 
the far spanwise boundary perpendicularly. Various combinations of spanwise grid point 
spacing have been studied. Both uniform and cosine-distributions on the lifting surface and 
uniform and stretched distributions outboard of the wingtip have been used. Figure 54 indicates 
the grid distribution for the case studied in ref. 120, the RAE tailplane model which was tested 
for pitching oscillations. It has a leading edge sweep angle of 50.2 deg., a taper ratio of 0.27 and 
a 10 percent thick NACA 64A010 section profile. Figure 55 shows comparisons of steady 
pressures for M = 0.80. The agreement is reasonable for the inboard stations. However, the 
agreement deteriorates at the outboard stations and the pressure expansion over the forward 
portion of the wing is generally underpredicted. The unsteady pressures for this case with k = 
0.490, fig. 56, _ show good agreement between calculated and measured trends and again 
deteriorating agreement outboard. 

The grids used for the calculations shown in figs. 53-56 contained 60 x 20 x 40 points in 
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the x, y and z directions for a total of 48,000 points. This is considered to be a medium-to- 
coarse grid and 5-6 times this number of grid points appear to be required to insure accurate 
calculations for a single lifting surface. In terms of computational efficiency, ref. 29 indicates 
a machine speed of 83 *is. per grid point per time step for a CDC 7600 computer. 

Sh ock-Generated Entropy . Gibbons et al.122 modified the XTRAN3S code to account for 
shock-generated entropy in a manner similar to ref. 119. The streamwise flux, eq. (31), is 
replaced by 


fi = (Y + 1)M 2 R (VV* - V 2 /2) + G<*> 2 ( 8 2 } 

and the pressure and wake condition modified as in eqs. 72-77. Reference 122 indicates that 
for a 3-D rectangular wing multiple solutions occur for aspect ratios greater than 24. For 
more reasonable aspect ratios, unacceptable solutions can be calculated if entropy generation is 
not taken into account. Figure 57 shows such a case wherein the unmodified code predicts a lift- 
curve slope twice as large as the modified code which is in agreement with an independent Euler 
equation calculation. Reference 122 also discusses implementation of Engquist-Osher monotone 
differencing in the XTRAN3S code. This resulted in improvement in calculated pressures near 
shocks. Figures 58 and 59 show original and modified XTRAN3S steady and unsteady pressures 
for the RAE tailplane model at M = 0.90 and k = 0.44. The entropy modifications result in 5- 
10 percent chord forward shifts of the shock and shock pulse in the outboard wing region. 

Treatment of Realistic Configurations 
All results presented thus far have dealt with unsteady aerodynamics for isolated lifting 
surfaces. In order to realize necessary improvements over existing aeroelastic analysis 
methods, computational methods are required to provide reliable predictions for complex 
configurations. 


Mu ltiple Lifting Surface Interference - Batinai23 extended the ADI algorithm of the 
XTRAN3S code to allow two lifting surfaces. Figure 60 illustrates a case of canard-wing 
interference in which unsteady loads are induced on the wing by the oscillating canard. This 

effect is obviously a function of the separation distance between the surfaces, the Mach number 
and the frequency. 
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Wing-Fuselage Interference - Batina 1 24,85 has also implemented modifications to the TSD 
equation to treat wing-fuselage interference. For a body at angle of attack ab and at yaw angle 

Pb, the exact flow tangency boundary condition may be written as 1 24 

N , + N „ ( 1+ V + N y<VPb> + N z<f* + “b> = 0 (83) 

where N(x, y, z, t) = 0 defines the body surface. Computationally, bodies are modeled by 
applying simplified boundary conditons on a prismatic surface rather than on the true surface 
as shown in fig. 61. This method is consistent with the small-disturbance approximation and 
treats bodies with sufficient accuracy to obtain the correct global effect on the flow field without 
the use of special grids or complicated coordinate transformations. The appropriate boundary 
conditons imposed on the computational surface are 

Upstream face: <|> x = - 1 ( 8 4 a ) 

Downstream face: <|> x = V exit - 1 ' (84b) 

♦, = c .tr- + iri- c A (84c) 

y y 
N N 

(84d) 

z z 


Left/right faces: 


Top/bottom faces: 


where Viniet and Vexit are inlet and exit flow velocities, respectively, specified in the case of a 

nacelle. The parameters Ct and Ca are thickness and angle-of-attack corrections derived from 

slender body theory to account for the differences between the true and computational body 
surfaces. 

Figure 62^4 gives results from these modifications to the XTRAN3S code for the RAE wing- 
fuselage tested at M = 0.90 and a = 1°. Results for the wing-alone case, not shown, indicate 
that fuselage effects are similar to the results shown for M = 0.91. Figure 63 shows the 
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calculated interference effect for this configuration for an assumed wing bending mode. The 
interference effect on the integrated generalized force, important for aeroelastic analysis is 
approximately 5 percent of the total. 

Approximate-Factorization of the TSD Equation 
The numerical stability restrictions of the ADI solution method for the 3-D TSD equation 
(29) have limited its applications for transonic unsteady calculations. An alternative solution 

algorithm 84 , termed an approximate-factorization (AF) method, has shown much improved 
stability characteristics. This AF method involves a local time linearization procedure coupled 
with a Newton iteration technique which is based upon the work of Shankar et al.t25 who applied 
the method to the full potential equation. It is formulated by first approximating the time 
derivative terms (<j>tt and <>xt) in eq. (30) by second-order accurate formulae, followed by the 

substitution of (|> = <}>* + A<> and the neglect of squares of derivatives of A<j>. In this method, <t>* is 

the current estimate of <|> n+1 which will be converged to the true potential <j> n+1 thus driving a<|> 
to zero. Performing these operations and summing terms results in 

2L A* + — A* - JL (E A6 + 2F $ A* + 2G<>* A<j> ) 

At 2 2 At ox * x x y y 


.n-1 ,n-2 


• 57 w >, + H < + H < ^ = - A 2 » - ^ - * 


At 


3<!)* - 4<t> n + <t> n l 
B x _JL " x 

2At 



(85) 


The right hand side is simply the TSD equation which may be evaluated using known potentials 

# 

<|>\ 4> n , . and 4> n -2. Transforming to computational coordinates, rewriting eq. (85) in 

conservation form and approximately factoring the left hand side into a triple product of 
operators yields 
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L L L A<{> = -R (<!>*, <T\ f' 2 ) 

s ^ s 


where 


( 86 ) 



, 3B K A 3 „ At 3 

+ 4A ^ X ^ 3£, 


F l— 


L 

T1 


^Li F ± 

x ^ 3ri 2 3r] 


(87a) 


(87b) 



^LAf± 

’* ^ 3C 3 3C 


(87c) 





(87d) 


F 2 *f (1 + H ^V 

<U 4 


(87e) 


F 3 ~T ( 8 7 f ) 

. At 2 , A 2<J>* - 5<J>" + V' 1 - «> n ' 2 3<t>: - 4<()” + (J)"' 1 

R = '^2A 2 * B -^ ^ L " 

^ At 2At 



[E^ + Fq 2 (J)* 2 + G(^ y 4>* + 0* ) 2 + 5l + ((,*) + (£, y 4>* + <j>* )] 
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+ — [— (qA + <t>* ) + H(j)* £ v 4>t + <!>*)] + [ 7 “ W ) ( 8 7 g ) 

an y * n ^ y s 11 3C ^ 


Equation (86) is in the form of Newton’s method for the solution of the nonlinear system of 
equations 


R (<)) n+1 ) = 0 


( 88 ) 


which is given by 


R (<t> n+1 ) = R (40 + (— ) . (4> n+1 - 40 = 0 

a<() 


(89) 


or 



* 3R -i * 
= 4> -[(— ) .1 R(<t>) 

a<i) *=* 


(90) 


By iteratively solving (90), A<*> will approach zero so that <|> n+1 = <)>*. In general, only one or 
two iterations are required to achieve acceptable accuracy since the Newton iteration process is 
quadratically convergent. Equation (86) is solved using three sweeps through the grid by 
sequentially applying the operators L4, L^, and as 

sweep: L„ A<j) = - R (91a) 


L A<J) = Ad> 
ri 


t| - sweep: 


(91b) 
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C - sweep: 


A$ = A<J> 


(91c) 


Central difference formulas are employed for all of the derivatives on the left-hand side of 
eqs. (91) except for the second term in the L4 operator (from the <)>xt term) which is backward 

differenced to maintain stability and the third term in the L= operator which is split into 

streamwise and spanwise components. The resulting terms are centrally differenced at subsonic 
points and the streamwise terms are upwind-biased at supersonic points using the Murman- 
Cole type dependent mixed difference operator. The terms on the right hand side of the £-sweep 
are also approximated using central-difference operators except for the <|>xt term which is 
backward differenced and for terms in the streamwise direction which are upwind biased at 
supersonic points. Since the L4, L^, and U; operators only contain derivatives in their 

respective coordinate directions, all three sweeps may be vectorized. Finally, ref. 126 
describes two additional improvements to the algorithm: second-order accurate supersonic 
differencing for improved accuracy and Engquist-Osher monotone differencing which again 
enhances stability. 

Figures 64 and 65, from ref. 125, compare results from the XTRAN3S ADI code and the AF 
algorithm for the F-5 wing model. Figure 64 compares steady pressures at M = 0.9 and a = 0° 
while fig. 65 compares unsteady upper surface pressures for k = 0.137. Results from the two 
algorithms agree, as they should, and they also are in very good agreement with experiment. 

The AF algorithm is also capable of calculating supersonic freestream conditons. Figure 66 
gives an example of unsteady lower surface pressures for the F-5 model at M = 1.1 showing 
good comparison with experiment for this relatively low supersonic case. 

The AF algorithm is implemented in a computer code termed CAP-TSD (Computational 
Aeroelasticity Program - Transonic Small Disturbance) developed at NASA Langley Research 

Center84. The code permits the aeroelastic analysis of complete aircraft through the modeling 
of multiple lifting surfaces and bodies including canard, wing, tail, control surfaces, launchers, 
pylons, fuselage, stores and nacelles. Reference 84 presents results for five configurations 
illustrating this capability. Figure 67 indicates the modeling of a General Dynamics F-16C 
aircraft model using four lifting surface and two bodies. Figure 68 compares calculated and 
measured steady pressures for this model at M = 0.9 and a = 2.38°. The agreement is 
considered very good. In fig. 69, calculated unsteady pressures for the.complete aircraft and for 
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the wing-alone are compared for the case of rigid pitching oscillations for k * 0.1. There is a 
general lack of unsteady experimental data for complex configurations with which to validate 
such computations. The grid used for these calculations contained 324,000 points. The 
calculations required 0.88 CPU seconds per time step or 2.7 ps. per grid point per time step on 
the CDC VPS-32 computer. Thirteen million words of memory was required. 

A primary application of unsteady aerodynamics for lifting surfaces is the calculation of 
aeroelastic response; that is, the response of elastic structures interacting with the 
aerodynamic loading. The integrity of the structure must be insured under all possible 
operating conditions and thus the prediction of instability, or flutter, boundaries is required. 
Reference 127 describes such calculations for a proposed AGARD standard aeroelastic flutter 
model configuration. This 45° sweptback wing with a taper ratio of 0.66 is shown in fig. 70. A 
finite element structural dynamic model provides normal vibration mode shapes, frequencies 
and generalized masses required for a flutter analysis. These uncoupled modal equations of 
motion are implemented in the CAP-TSD code and time-accurate transient response calculations 

obtained 127 . 

Figure 71 shows the mode shapes of the four lowest frequency modes used in the analysis. 
Transient responses for values of dynamic pressure below and above the flutter boundary are 
processed to obtain the frequency and damping of the aerodynamically coupled modal responses. 
Figure 72 compares experimental and computed flutter boundaries for Mach numbers from 
0.338 to 1.141. Figure 72(a) shows the flutter speed index, U/bco a Vp, and fig. 72(b) the 

frequency ratio, of the flutter mode, where wa is the first torsion mode wind-off 

frequency. For this 4 percent thick wing, transonic effects are delayed to high subsonic Mach 
numbers and linear theory results from both CAP-TSD and a kernel function program, are in 
very good agreement with experiment up to M = 0.98. The three nonlinear CAP-TSD subsonic 
flutter calculations better agree with experiment than the linear theory, particularly for the 
change in slope of the flutter boundary near M = 0.95. Note the excellent prediction of the 
supersonic "backside” of the flutter dip. 

Unsteady Potential Equation 

The length of the previous sections dealing with solution methods and modifications for the 
unsteady TSP potential equation reflects the large amount of effort which has been expended 
upon this equation level. Fewer results are available for the higher equation levels. However, 
as work on the TSD equation has obviously matured, an acceleration of work on the potential 
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equation, and the Euler and Navier-Stokes equations is evident. In this and the remaining 
sections, key algorithms and results for these fluid dynamic models will be summarized. The 
major issues to be addressed involve computational grids and increased problem size. All of the 
fluid dynamic flow models above the TSD equation require solutions to be obtained upon the 
actual body surface, which is usually accomplished with a body conforming grid. For problems 
involving body motion, such grids probably need to be moved aiong with the body. 

The 2-D unsteady potential equation, eq. (23), is solved in ref. 82 using an approximate 
factorization algorithm. The time derivative of density is linearized about previous time levels 
such that conservation form is maintained and the resulting equation becomes 


where 


.2 

L = [I + hU n 8 8 (p Aj)" 8 ] 

% s p n s s 

.2 

L = [I + hW n 8 - — — 8 (p A 3 ) n 8 ] 

; ^ p n ; ; 

F = (O n - d>" *) + -fi— (O n - 2<t> n 1 + O n 2 ) + — (p - P ) 

+ h (U n l 5 5 + W"- 1 5 ; ) (4>° - *"■') + ^ [S t (^■) n + (pw>“ - FJ 

The terms Ai and A 3 are metric terms, 6<; and 5^ are central difference operators and 
AO = o" 1 - O", h = At, p-p/J, P = P 2 ' 7 « 

Reference 82 uses the concept of flux-biasing'*® to stabilize Ihe numerical scheme. The 

C ' *-( 
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spatial terms on the right hand side of eq. (92) are centrally differenced about node (i, j) to 
give 


rPU x _,PU 


5 (^-). . = (i2i) . 

4 J '«.J V J M+l/2,j ' j •'i-l/2,j 


(93a) 


S ? (pW) y = (pW). J+1/2 -(pW).j. 


1/2 


(93b) 


The biased density in the 4 direction is defined as 


i+!/2.j q i+J/2 j *- Pq ' ^ pq ^i+i/ 2 .j + (P^i-l/i.p 


where 


(pq)* = pq - 


P q 


q > q* 


o 


q<q 


This flux biasing has the effect of introducing artificial viscosity in supersonic regions which is 
necessary to capture shocks. Although it is generally necessary to bias the density in both 
computational directions, ref. 82 found that biasing only the 4 direction was satisfactory. 
Reference 128 explores the connections between this flux biasing and the Murman-Cole and 
Engquist-Osher differencing schemes. 

The flux-biasing solution method has the following desirable features: (a) it accurately 
tracks sonic conditons and requires no empirical constants to specify the amount of artificial 
viscosity, (b) it produces no velocity overshoots at shock waves, allowing for larger time steps 
for unsteady calculations, (c) it produces well-defined, monotone shock profiles with a 
maximum two point transition between upstream and downstream states, and (d) it dissipates 
expansion shock waves, ruling out solutions with such nonphysical characteristics. Shankar et 

al. 130 have also used flux-biased differencing in unsteady potential equation calculations. 
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Reference 82 also presents an entropy correction method to account for shock generated 
entropy. The correction consists of replacing the isentropic density by a nonisentropic density 


p = p.e 


-As/R 


(34) 


where the entropy jump as is a function of the Mach number normal to the shock. The 
nonisentropic pressure coefficient is given by 


C„ = 

y 


2 -As/R. y - , 

— — [(p ; e )'-l] 



and, as for the TSD equation nonisentropic modifications, the entropy jump is convected along 
the wake downstream of the trailing edge. 

Treatment of moving, body conforming grids has been handled in several different manners. 
For 2-D airfoil sections and 3-D rigid wings, a single grid generated for a nominal body 

orientation may be calculated. Rotations and translations of the entire grid 81 can then be used to 
track the motion of such bodies. An alternative is to calculate two grids for the extremes of body 

motion and linearly interpolate grid point locations for intermediate body orientations 26 . 

Figure 73 82 illustrates this latter method showing calculated grids for 0° and 45° and 
interpolated grids for 15 deg. and 30 deg. The large amplitudes of this examples serve to 
demonstrate the method. The double wake grid line for these potential equation applicatons is 
shown opened for clarity. For the potential equation, the location of this wake cut is important 
since it defines the path taken by the convected vorticity. Reference 131 studies the effect of 
using linear, quadratic and cubic curves to define the wake cut and shows a singificant effect 
upon calculated lift results. 

Figure 74, from ref. 82, shows isentropic and nonisentropic potential equation calculations 
for the NACA 0012 airfoil oscillating in pitch at M = 0.755 and for a(x) = 0.016° + 2.51° sin 
(kx) where k = 0.814. The effects of the entropy corrections are to weaken the shock and move 
it forward, in bettter agreement with the experimental data. At points in the cycle where the 
shocks become strong, the measured pressures immediately behind the shocks show the effects 
of boundary layer thickening, which is not included in these inviscid calculations. Reference 
129 gives further examples of potential equation calculations for the AGARD Standard 
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Configurations. 


Euler and Navier-Stokes Equations 

Since the Euler equations, (14), may be obtained from the Navier-Stokes equations by 
deleting viscous terms such as Hv in eq. (1), solution algorithms for both equations sets may be 

discussed together. Edwards and Thomas 1 survey methods which have been used. 

The time-dependent Euler equations form a hyperbolic system of equations, and much of the 

recent progress in algorithm development 132 ’ 137 has hinged upon the incorporation of the 
signal propagation features of the differential equation into the numerical algorithm. There are 
several methods of incorporating this information into a difference scheme, for example flux- 
vector-splitting or flux- difference-splitting, and an excellent review of the current 
developments in the field is given by Roe in Ref. 138. The advantages of incorporating an 
upwind-biased discretization into a numerical algorithm are twofold: (1) the scheme becomes 
naturally dissipative so that no adjustable constants need to be fine-tuned to a particular 
application and (2) improved implicit schemes can be devised for more efficient solution to 
both steady and time-dependent problems. Both of these advantages offset the disadvantage that 
approximately twice as many operations per time step are required to implement an upwind 
scheme as opposed to a central difference scheme. 

Most of the calculations made to date with upwind difference schemes, especially for 
airfoils/wings, have been steady-state applications, for which comparable accuracy can be 
obtained by central difference methods with added artificial viscosity. The advantages of upwind 
differencing should be more significant for time-dependent problems, however, where the 

ability to treat rapid movement of flows with shocks is required. Roe 138 gives several 
examples of shock-propagation computations in two-dimensions which demonstrate clearly the 
advantages of a characteristic-based scheme. Viscous effects can also be readily introduced into 
upwind difference schemes developed for the Euler equations by central differencing the shear 

stress/heat transfer terms 139 * 140 . 

The time-accurate computations made by Steger and Bailey 25 and Chyu et al. 26 - 27 used a 
spatially-split approximate-factorization (ADI) scheme, which is unconditionally stable in two 
dimensions but at most conditionally stable in three dimensions. Alternate factorizations are 
possible with the incorporation of an upwind difference discretization in one or more coordinate 

directions which can lead to unconditionally stable 3-D algorithms 132 . A two-factor eigenvalue 
split scheme for the Euler equations has an increased stability limit and fewer operations than 
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the spatially-split scheme, although the operations are not completely vectorizable. Belk 4 2 
computed steady and time-dependent inviscid flows for the NASA RSW model with such an 
algorithm in combination with a blocked-grid strategy. Ying et al. 14 2 used upwind differences 
in a single coordinate direction and constructed a two-factor unconditionally stable algorithm 
for which thin-layer viscous effects are readily incorporated. Applications of the thin-layer 
Navier-Stokes equations to the high-angle-of-attack unsteady flow over a hemisphere-cylinder 

are made 142 . Several of ihese alternate factorizations are investigated in ihe context of 

efficient algorithms for three-dimensional steady-state problems by Anderson et al. 14 3. 144 . 

The use of multigrid techniques to accelerate convergence to the steady-state is becoming 
widespread in the aerodynamic community. These techniques can also be used for time- 
dependent flows. For instance, multigrid techniques could be used to efficiently solve the large 
banded matrix equations arising from implicit time discretizations, the solution of which is 

generally approximated through relaxation and/or factorization methods. Jesperson 14 5 has 
demonstrated a time-accurate multiple grid procedure which was used to overcome the small 
timestep limitation of an explicit scheme. With the growing memory of today’s computers (the 
Numerical Aerodynamic Simulator has 256 million words of memory) it becomes feasible to 
solve the banded matrices by direct Gaussian elimination, rather than by approximate 
techniques. The structure of future implicit algorithms for both steady and time-dependent 
problems will likely involve a multiple grid algorithm with direct elimination techniques used 
on the coarser grid levels. 

Anderson et al.81 implement the solution to Euler equations for a moving grid, eq. (14), 
using flux vector splitting with upwind differencing. An iterative Newton linearization is used 
to advance the solution in time similarly to eqs. (88) - (89). That is, eq. (1) is reformulated 
as 


L(Q" +1 )=0 (96) 

where the form of the operator is 

[(1 +4)Q n+1 -(l+<l))Q n + 4Q n - 1 ]-L 

2 2 At 
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(97) 


If ()> = 0, the scheme is first order accurate in time while if <i> = 1 . the scheme is temporally 
second order accurate. Conditions required for the split flux terms in eq. (97) are given in ref. 
81. Equation (96) is a nonlinear equation which can be solved iteratively by Newton 
linearization 


3L 


u+i 


(Q - Q) = - L(Q) 


(98) 


A A A A 

where I is a sequence of iterates and at convergence, Q l+1 - Q ,=0 , Q 1 = Q n+1 . A spatially-split 
approximate factorization scheme is used to solve (90) in three sweeps in each coordinate 
direction 


[(I + $) + At8 A* + AtS^A ]AQ* = - AtL(Q l ) 
2 \ k 


[(I + £) + At5B + + AtS^'lAQ** = (I + ^)AQ* 
2 n 'n 2 

[(I + 1) + AtS^cf + At5^]AQ 1 = (I + |)AQ** 



A 1 A 1 
= Q + aq 1 


AAA AAA 

where A±, B±, C± correspond to Jacobean matrices of F±, G±, H± respectively. 

Figures 75 and 76 from ref. 81, show calculations for the pitching NACA0012 airfoil for 
the same conditions as in figures 37(d) and (73). The computations were obtained on a 193 x 
33 C-grid using a time step of 0.10. Figure 75 compares results from two alternative splitting 
methods; flux vector splitting, FVS, and flux difference splitting, FDS. The comparison of both 
of the euler equation results with experiment is very good except near the base of strong shocks. 

Reference 146 presents the Euler equations in a formulation wherein the fluid dynamic 
equations are developed for a coordinate system rigidly attached to the translating and rotating 
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body. Solutions for two-dimensional oscillating airfoil problems are obtained using an implicit 
approximate factorization method with artificial dissipation. Results for the oscillating 
NACA0012 airfoil are shown in figs. 77 and 78. A C-grid with 128 x 64 cells was used for 
these calculations. Comparing figs. 38(d), 76, and 77 indicates that the two Euler equation 
results appear to be converged while the TSD solution indicates the magnitude to be expected 
from viscous effects. 

Reference 146 also presents results for unsteady locally conical flow about a sharp-edged 
delta wing in supersonic flow for rolling oscillation about zero angle of attack. Figure 79 shows 
the spanwise pressure distribution at four instants during a cycle of oscillation for M = 2 and k 
= coc/U = 1.337. The formstion of the Seeding edge vortex end its migration ere o!eer!y evident 
as is the phase lag of the pressure loading. 

The final example of Euler equation calculations to be presented is from the flux-vector 
splitting scheme of ref. 81. Figures 80 and 81 show steady and unsteady pressures for the F-5 
wing model oscillating in pitch. Figure 80 presents comparisons with experimental steady 
pressures for Mach numbers from 0.90 to 1.328. Agreement is very good to excellent except 
near the strong shock at M = 0.95. Unsteady pressure comparisons, fig. 81 , for M = 0.95 and 
1 .32 show very good agreement in pressure levels. There are generally insufficient data to 
resolve the shock pulses near the leading and trailing edges at M = 0.95. These calculations 
were obtained on a 129 x 33 x 33 C-H mesh. A time step of 0.05 was used requiring 
approximately 240 time steps per cycle of oscillation. 

The retention of viscous terms leads back to the Navier-Stokes equations represented by eqs. 
(1) - (13). The detailed viscous flow modeling of which this equation set is capable makes it 
appropriate for the study of viscous dominated unsteady flows characterized in fig. 3 as type II 
and type III. They may aso be used in the calibration of lower equation level flow models for 
appropriate classes of unsteady fows such as type I attached flow. 

Steger and Bailey 25 provided an early computational demonstration of the use of CFD 
methods in aeroelasticity. They studied a case of aileron buzz for the P-80 aircraft which had 
been tested in a wind tunnel. Aileron buzz is a one-degree-of-freedom aeroelastic instability, 
usually of limited amplitude, which may be encountered for transonic flow conditions. They 
implemented the Beam-Warming implicit Approximate Factorization (AF) solution algorithm, 
using an algebraic eddy viscosity turbulence model. A novel treatment of the computational grid 
was used to follow the aileron motion with a conforming grid. A simple shearing transformation 
in the coordinate normal to the aileron was used. 
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Figure 982shows the limit amplitude "aileron buzz" oscillation which was calculated for M 
= 0.82, Re = 20 x 10 6 and a = -1. The calculation was initiated with a 4 degree aileron offset. 
This and other calculations successfully reproduced the experimental buzz boundary. The 
computed reduced frequency was k = 0.38. The shock motion observed in the calculations was 
type B, and type II intermittent fow separation is shown in ref. 25. The code was capable of 
being run in an inviscid mode (EE mode) and several such calculations were made. Below M = 
0.84 the aileron exhibited damped oscillations of about k = 0.36 whereas divergent oscillations 
(k = 0.39) were calculated at M = 0.84. Hence the tendency to oscillate at a given frequency 
derives from the inviscid flow equations while the viscosity apparently plays the key role of 
limiting the amplitude of oscillation. These calculations were performed on a 76 x 42 grid and 
required approximately 1.5 sec of CDC 7600 computer time per time step or 460 ps. per grid 
point per iteration. Nondimensional time steps of 0.005 - 0.01 were used (based on chord). 

Chyu and his coworkers 26 ’ 2 ^ used this same method along with the grid interpolation 
method described above to study the moderate shock case 26 and the shock-induced separation 
case 27 for pitching oscillations of the NACA 64A010A airfoil shown in figs. 83 and 84. 
Comparison of the interactive viscous-TSD (IV-TSD) results of fig. 35 with the Navier-Stokes 
results of fig. 83 is instructive. The two sets of calculations are in very good agreement for this 
type I flow condition. As for the type II - III flow condition shown in fig. 84, note that the full 
and thin-layer Navier-Stokes results show no significant differences except near the shock 
where the difference is not large. Thus the thin-layer Navier-Stokes equation appear to be 
viable for this class of unsteady flow problem. These calculations were obtained on a 139 x 49 
C-grid using 2620 steps per cycle of oscillation and the time per step on the CRAY X-MP 
computer was: full NS, 0.33 sec; TL-NS. 0.17 sec; EE. 0.17 sec, corresponding to 25 - 44 
psec. per grid point per time step. 

Rumsey and Anderson78 describe an extension of the flux-vector split, approximate 
factorization upwind scheme described in eqs. (96) - (98) to the 2-D thin-layer Navier- 
Stokes equations. This method is developed for body conforming moving grid systems and is 
given in eqs. (1) - (13). Figure 85 gives their solution for the oscillating NACA 0012 airfoil 
at M = 0.6 and k = 0.081. This case has already been presented for IV-TSD calculations in fig. 
38(c) and for EE calculations in fig. 78. Figure 85 shows the effect of grid refinement 
indicating little change in going from a 129 x 29 grid to a 257 x 97 grid. The EE and TL-NS 
solutions appear to be quite similar while the TSD solution shows some differences for this high 
lift case. This computational algorithm is highly vectorized for use on either the CDC CYBER 


63 


205 or the CRAY 2 computers. Computational speed averages approximately 16 psec per grid 
point per time step and the memory required in kilowords is about 0.260 x (meshsize). 


Periodic Aerodynamic Oscillations 

In the previous sections, algorithms and experimental results directed at moving or 
oscillating lifting surfaces have been presented. There is an important class of experimentally 
observed unsteady flows wherein periodic flow is observed over very narrow ranges of test 
conditions for perfectly rigid bodies. These flows tend to be found in transition regions between 
attached and separated flow conditons and are recognized as highly sensitive data which may be 


used for the validation of computational methods 1 . 

In order to provide experimental data for validation of viscous flow CFD computer codes, 

McDevitt et al. 43 conducted tests on a rigid 18 percent thick circular arc airfoil. Figure 86 
illustrates the parameters of the experiment which was designed to encounter both trailing- 
edge and shock-induced separations at high Reynolds numbers within the wind tunnel operating 
limits. Over a narrow range of Mach number, 0.73 < M < 0.78, oscillatory flow separation was 


observed, fig. 87 43 . McDevitt 44 states that the oscillations involve predominantly type C shock 
motion with small regions of type A motion near the onset of the periodic oscillations, fig. 88. 
The reduced frequency of the oscillations is k * 0.48 for a = 0° and varies little with angle-of- 
attack. 


Levy 4 5 successfully computed such oscillations for this airfol using a Navier-Stokes flow 
solver. Levy’s code was a. modification of the code of ref. 148 and used MacCormack's explicit 
solution scheme with an algebraic eddy viscosity model. Levy modified the code to simulate the 
contoured wind tunnel walls. Figure 89 shows steady computed Mach contours for Mach 
numbers of 0.72 and 0.78, corresponding to trailing-edge and shock-induced separations, 
respectively, and unsteady flow with oscillatory trailing-edge/shock-induced separation for M 
= 0.754. The reduced frequency of the computed oscillations is 0.40, about 20 percent lower 
than the measured frequency. Note particularly the lower surface Mach contours of the third 
frame for M = 0.754. The few lines indicate the collapse of the supersonic region for this 
portion of the cycle. Also note the dimpled nature of these Mach lines under the airfoil surface. 
These features will be discussed in more detail below. 

Subsequent tests on circular arc airfoils of thicknesses from 10 to 20 percent were 

performed by Mabey 46 , obtaining similar periodic oscillations. The Mach number range of the 
oscillations increases with decreasing thickness as does the oscillation frequency, remaining in 
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the range of 0.4 < k < 0.55 depending on thickness and wind tunnel wall condition. Mabey et 

al. 4 ^ and Levy48 gj ve detailed comparisons of Navier-Stokes calculations with experiment for a 
14 percent thick airfoil, reproducing qualitatively the details of the oscillatory flow. 

These unsteady periodic flows encountered over limited ranges of Mach number and 
triggered by oscillating trailing-edge/shock-induced boundary layer separation are just 
recently coming within the range of computational methods. The weakest link for this capability 
is the uncertainty in the turbulence modeling of complex separated flows, since rapid progress 
continues to be made in the development of improved algorithms and faster computers. 

It is well known that separated flows depart strongly from equilibrium - type behavior, so 
that at a minimum some account of the non-equilibrium "upstream history" effects should be 
included in the computations. Some encouraging results along this line have been obtained by 
LeBalleur149 with an integral boundary layer model and JohnsonlSO with an eddy-viscosity 
Reynolds-shear stress closure model. Simpsonisi recently reviewed calculation methods for 
turbulent separated flows and C‘oakleyi52 compared several methods for airfoil applications. 

On the other hand, Levy45 wa s able to reproduce the unsteady periodic flow behavior of the 
1 8 percent circular-arc airfoil using an equilibrium two-layer algebraic model. The steady 
flow at Mach numbers below the range of periodic flow, characterized by trailing-edge 
separation, was predicted accurately. Levy demonstrated that the influence of the channel walls 
had a substantial impact on the comparisons with experiment, especially at Mach numbers away 

from the design point. This effect was not considered in the earlier comparisons of Deiwert148 
with the experimental results. The steady flow at a Mach number above the range of periodic 
flow, characterized by shock-induced separation, was not accurately predicted, as the 
calculation demonstrated a normal shock pattern (fig. 89) with trailing-edge pressure 

recovery, whereas the experiment indicated an oblique shock pattern and a constant pressure 
region downstream of the shock. 

In addition to Levy's calculations, the unsteady periodic behavior for the 18-percent 
biconvex airfoil has also been computed by Stegeri53 a nd by LeBalleurl49. steger's calculation 
was for an airfoil in free-air with an implicit Navier-Stokes code using the Baldwin-Lomax 
algebraic model. The unsteady flow occurred at a higher Mach number (M = 0.783) than that of 
Levy (M = 0.754), which can partly be attributed to the free-air boundary conditions. The 
computed reduced frequency (0.41) was remarkably close to that of Levy (0.40) although both 
are low in comparison to experiment (0.48). LeBalleur's recent calculations were also made in 
free-air with a small disturbance potential method including an interacted two-equation 
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integral viscous model. Steady shock-induced separation was computed at M = 0.788 and 
unsteady periodic flow at M = 0.76. The reduced frequency (0.34) was lower than either of the 
two Navier-Stokes solutions. 

Some calculations have been made using the implicit upwind-biased Navier-Stokes 
algorithm described in ref. 154 using an algebraic turbulence model. The tunnel walls were 
modeled and boundary conditons appropriate for internal flow were used, i.e., the downstream 
pressure and upstream enthalpy, entropy, and flow direction were specified. The results 
indicated unsteady flow at a higher Mach number than Levy; steady trailing-edge separation 
occurred at M = 0.754 and unsteady periodic flow at M = 0.78, although the Mach number for 
Anpftt fho imctoaHlnocc wac concitivo tn whpthpr nr nnt thp divprnpnrp of tho tunnol boundary 
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to account for boundary layer growth was included. Figure 90 shows Mach contours through one 
half-cycle of oscillation (near maximum lift to minimum lift) indicating the forward 
movement, disappearance, and subsequent formation near the trailing edge of the lower surface 
shock. The reduced frequency of the type B unsteady motion was 0.406, in close agreement with 
the calculations of both Levy (compare fig. 89) and Steger. The implicit calculations were made 
with a time step of 0.01 and a computational time of 18 ps. per grid point per time step on the 
CYBER 205 computer. 

The calculation of the unsteady periodic flow boundaries for airfoils is a fruitful area for the 
development and validation of computational methods. Experimental pressure data 4 3-47 over a 
wide range of Reynolds number is available, although detailed boundary-layer measurements 
are not. For the 18-percent biconvex experiments of McDevitt, a substantial hysteresis effect 
in the unsteady flow boundary was found. This aspect has not been demonstrated with 
computational methods as yet, but it would be expected, based on the above discussion, that 
computational modeling as close as possible to that of the experimental conditions will be a 
critical consideration. The most interesting behavior, and the most challenging from the 
computational viewpoint, occurs in the transitional region from laminar to turbulent flow. In 
the experiment of McDevitt 4 3, the Mach number range for the observed unsteady flow 
diminished near a Reynolds number of 3 x 10® (fig. 87) and in the experiments of Mabey 47 , it 

disappeared completely in the range of Reynolds number from about 3 x 20 6 to 5 x 10 6 . Scale 
effects in either experiment were not significant once turbulent flow is fully established ahead 
of the shock. 

The frequency of these oscillations is of interest in that the flow mechanism causing the 
unsteadiness might be identified via the characteristic time constants of signal propagation 
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within the various flow regions. Tijdeman3 noted an almost linear relation between the phase 
lag of the shock motion and the airfoil motion for type A shock motion with a well-developed 
shock (for pitch oscillations of the NLR 7301 airfoil). He related this to the signal propagation 

time from the trailing-edge to the shock. Mabey 47 , commenting on characteristic time 
constants for the 14 percent circular arc airfoil periodic oscillations, notes that this reasoning 
leads to reduced frequency parameters of 1.15 to 1.8, much higher than the observed 
frequencies. 

Three items mentioned above germane to this discussion are; Steger and Bailey's 25 inviscid 
EE calculations for the aileron buzz case, LeBalleur and Gerodrous-Lavigne's 149 interacted 

viscous-TSD code result for periodic oscillations, and Batina’s 108 demonstration of the 
possibility of aerodynamic resonance with an inviscid TSD code. The occurrence of damped and 
undamped oscillations observed for inviscid flows at nearly the same frequency as the 

oscillations in viscous flow 25 (k = 0.36-0.39) implies that the flow mechanism determining 
the oscillation frequency derives from the dynamics of the inviscid flow region. Furthermore, 
the results of Refs. 108 and 149 give impetus to studying this effect with a TSD code. 
Accordingly, calculations were made with the XTRAN2L code of the aerodynamic response for the 
18 percent thick circular arc airfoil due to trailing-edge 25 percent chord flap motions. The 
nonisentropic modifications described in eqs. (70) - (77) were used to obtain solutions for 
this strong-shock case. Figure 91 presents the resulting Cm5 frequency response function for 

Mach numbers of 0.66-0.74. There is a very marked development of an aerodynaic resonance 
effect as Mach number increases. For M = 0.74 the airfoil resonance frequency is k = 0.32, 
very cose to the computational conditions of ref. 149 for periodic oscillations (M = 0.76 and k 
= 0.36). 


CONCLUDING REMARKS 

These lectures have summarized the status of computational unsteady aerodynamic for 
lifting surfaces. The fluid dynamic flow models appropriate to the several levels of physical 
models available have been presented along with details of solution algorithms. The subject has 
been differentiated with respect to the difficulty of flow modeling and computational 
requirements by distinguishing unsteady flow types as: I, attached flow; II, mixed attached and 
separated flow; and III, separated flow. 

For type I flows, computational methods have matured with a steady progression of improved 
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techniques for flow simulation. Significant efforts have been devoted to understanding the 
effects of equational level, computational grid, boundary conditions, and interactive viscous 
modeling. Extensive comparisons with experimental data sets have been made with small 
disturbance potential (both linear and transonic), potential, Euler and Navier-Slokes equation 
solvers for two-dimensional airfoil cases and an understanding of the range of validity of the 
various methods can be made. Less extensive comparisons are available for wings and even 
fewer comparisons for complex configurations involving interference effects. Experimental 
unsteady data is needed for such configurations in order to validate computational methods which 
can now treat complete aircraft for transonic flow conditions. 


Progress in solution algorithms for unsteady aerodynamic has been significant with large 


decreases in required computer resources due to larger time steps allowed by more stable 
solution algorithms. The treatment of body conforming grids for deforming aeroelastic vehicles 
needs further attention to fully utilize the computational methods available. 

Computational aeroelastic analysis has demonstrated capability for prediction of complete 
transonic flutter boundaries for wings, including significant "transonic dip" features. Many 
more transonic flutter calculations will be needed to fully validate computational methods for 
transonic flutter predictions as critical features of these stability boundaries frequently 
involve difficult flow conditions, such as type II mixed flow. 
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Table 1. References giving comparisons of experimental and calculated 
two-dimensional unsteady aerodynamics 
(References from Edwards and Thomas - !) 

TSD FP EE NS 


NACA 64A006 

6. 103*. 88. 89. 87* 

102, 116 

5. 15 

1 

1 

NACA 64A010A 

80 . 82, 138 . 8C’. 89, 112*, 
87*. 93. 97+ , 114*. 105*-108* 

117, 116. 79, 119 

15, 30, 32 

19, ia ! 

i 

NACA 0012 

88, 90*. 109*, 97+ , 96+ , 93 

•119, 100+ 

29, 28, 31 

152, 44, lip 
43, 154 

NIR 7301 

89. 88. 107*. 114*, 97+ 

117 

16, 30, 32 


WB A-3 

112*, 89, 82, 138, 114*. 63 




Supercritical 

Airfoils 

154*. 109*, 65*. 66 




Circular arc 
Airfoils 

106*, 107*, 108* 



35, 37, 36, 40 

Other 



151 

17, 42, 43, 
148 


♦ Interacted boundary layer model 

♦ Noni sentropic corrections 


Table 2. References giving comparisons of experimental and calculated 
three-dimensional unsteady aerodynamics and aeroelasticity 
(References from Edwards and Thomas 1 ) 


TSO FP EE 


F -5 Model 

115. 157. 137. 23 

24, 136, 153. 102 


NORA 

76. 22 

64. 26 


• LANN 

100. 156 

45. 40. 24, 76 


RAE Wing A 


64. 43, 26. 136 


RSW 

116, 75 

40, 66, 76. 25 

96. 34 

RAE TAIIpUne 

IK, 85*. 155. 98* 


75. 33 

Other 

110*. 113*, 71*. 73* 

64*. 65*. 26*. 70* 



♦ Nonlsentroplc corrections 

* AeroeUstlc and flutter comparisons 





Table 4. Local singularity solutions for oscillating Table 5. Global singularities and situations where the local 

lifting surfaces. Note that in all cases the control flow may be irregular on interfering surfaces 

surface angle is 8(y,t) = 5(y)e*l“t (References given in Ref. 93) 

(Reference given in Ref. 93). 
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Integral, lief. 6) 













TABLE 6 

Current Levels of Accuracy for Aeroelastic Analysis 


Wing Flutter 

- 1 

0% 

Gust Response 

~ 1 

0% 

Loads 

5 - 

10% 

Control Effectiveness 

25 

- 50% 

Control Hinge Moments 

25 

- 50% 


Buffet Loads 


20 - 30% 



TABLE 7 

COMPUTER RESOURCE REQUIREMENTS FOR FLUTTER BOUNDARY 

(From Ref. 100) 


WING/BODY/CANARD CONFIGURATION 
10 MACH NUMBERS (40 CASES) PER ANALYSIS 

OPS OPS 

TIME = (GRID PTS) X (GRID PTS X ITER) X (ITER)/( SEC ) 


FLOW REGION 

FLOW MODEL 

VPS-32 
(100 MFLOPS) 

NAS 

(250 MFLOPS) 

1, MAYBE II 

TSD WITH 2-D 
STRIP BOUNDARY 
LAYER 

30 HOURS 

12 HOURS 

1, MAYBE II 

POTENTIAL WITH 2-D 
STRIP BOUNDARY LAYER 

45 HOURS 

18 HOURS 

1, II. MAYBE III 

EULER WITH 2-D 
STRIP BOUNDARY LAYER 

65 HOURS 

26 HOURS 


II, III 


NAVIER-STOKES 
(RE = 108) 


1611 HOURS 


644 HOURS 







TABLE 8 

COMPUTER RESOURCE REQUIREMENTS TO DETERMINE FLUTTER 
POINT AT A SPECIFIED MACH NUMBER 
(From Ref. 100) 

(4000 TIME STEPS PER FLUTTER POINT) 





CPU HOURS 

OONFK51 JRATiON 

FIOW MODEL 

ORID POINTS 

'VPS-32' 

COMPLETE AIRCRAFT 

TSD 

0.75M 

2.3* 

COMPLETE AIRCRAFT 

FULL NAVIER-ST OKES 
(RE = 10 MILLION) 

7.00M 

77.8** 


•BASED ON ACTUAL CASES 

••ASSUMES COMPUTATIONAL SPEED OF 100 MFLOPS 
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Equivalent 

airspeed 



Fig. 1 Graphical representation of minimum 
required flutter margin, (Ref. 2). 


M.. CONSTANT Cji ft • CONSTANT 



Fig. 2 Main characteristics and boundaries for 
the onset of separation for the NLR 7301 
supercritical airfoil. 

(Tijdeman, Ref. 3) 



Fig. 3 Characteristics of attached and 

separated flow for complete aircraft, 
(after Bradley, Ref. 5) 
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Fig. 4 Comparison of conventional wing flutter 
boundary with adjusted supercritical 
wing boundary 
(Farmer et al.. Ref. 6) 
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Fig. 5 DAST ARW-2 region of shock-induced 
high dynamic response, (Ref. 11). 



30K 


Altitude, 

feet 


0 



Mach number 


Fig. 6 High dynamic response region due to 
wing/store loading, (Ref. 100). 
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Angle at attack reduced toward 0* 


Uadi number 


High dynamic response region due to 
angle-of-attack, (Ref. 12). 


Fig. 9 AGARD Structures and Materials Panel 
two-dimensional standard aerceiastic 
configurations. 

(Bland, Ref. 16). 



Lift 

Coefficient 




Fig. 10 AGARD Structures and Materials Panel 
three-dimensional standard aeroelastic 
configurations. 

(Bland, Ref. 17). 


Fig. 8 Regions of vortex-induced buffet 
loads, (Ref. 100). 
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ig 13 Surface and panel geometry for 
doublet-lattice method. 

(Albano and Rodden. Ref. 95) 


RefGeneralized Force/q) 





Fig. 16 Pressure coefficient distribution due to 
airfoil plunging as a function of x and k 
at Mach 0.7. 

(Edwards, Ref. 98) 


Fig. 14 


Comparison of generalized force 
vectors predicted by RHOIV and 
doublet lattice on a rectangular wing 
of AR = 4, pitching about the trailing edge 
at k « 0.3, M = 0.9, wave number = 2.7. 
(Rowe and Cunningham, Ref. 96) 
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Fig. 15 Hinge moment vectors due to control surface 
motions for a wave number of 5.4. 

(Rowe and Cunningham, Ref. 96) 



Fig. 17 Pressure coefficient distribution due to 
airfoil plunging as a function of x and k 
at Mach 1.02 and k = 0.4 
(Edwards, Ref. 98) 
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Fig. 20 Upper surface pressure coefficients of an NACA 64A006 airfoil with oscillating 
trailing-edge flap showing type A, B. and C shock motions. 
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Fig. 23 Lift and moment coefficient frequency 
response functions from pulse transfer 
analysis. 

(Batina, Ref. 108) 



Fig.. 24 Unsteady forces for a 6% parabolic an 
airfoil calculated by pulse and 
oscillatory analyses; XTRAN2L default 
grid, M = 0.85. 

(Seidel et al„ Ref. 107) 


.05 



Fig. 25 Lift coefficient time histories due 
to pulsed pitching motion for flate 
plate airfoil at M = 0.85. 
(Edwards, Ref. 155) 
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RBC - reflecting boundary conditions 
NRBC - nonreflecting boundary conditions 




l 




REOUCED FREQUENCY. k 


a) RBC, 20c x 25c grid 


b) RBC, 200c x 2327c grid 



0 .5 1.0 15 r.o 


REDUCED FREQUENCY, k 

c) NRBC. 20c x 25c grid 

Fig. 26 Lift coefficient frequency response 
functions due to pulsed pitching 
motion for flap plate airfoil at 
M = 0.85. 

(Seidel et ah. Ref. 107) 
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Fig. 27 Steady pressure distribuitons for the NACA 64A006 airfoil 
(Howlett and Bland, Ref. 110) 
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Fig. 28 Effect of Mach number on boundary-layer 

displacement thickness for NACA 64A006 airfoil. 
(Howlett and Bland, Ref. 110) 
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Fig. 30 


Comparison of unsteady forces for the NACA 64A006 airfoil at 
(Howlett and Bland, Ref. 110) 

.05 


0 


So » 1° and k - 0.06 
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Fig. 31 Comparison of unsteady forces for the NACA 64A006 airfoil at So = 1° and k - 0.24 
(Howlett and Bland, Ref. 110). 
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(b) M = 0.796, Re = 12.5 x 106 . 



Ul M = 0.5 



Fig. 33 Effect of Reynolds number on boundary-layer 
displacement thickness for NACA 64A010A 
airfoil in steady flow. 

(Howlett and Bland, Ref. 110). 
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Fig. 32 Steady pressure distribution for the 
NACA 64A010A airfoil, am = 0°. 
(Howlett and Bland, Ref. 110). 


Fig. 34 Unsteady upper surface distribution for the 
NACA 64A010A airfoil at M = 0.5, ao = 1°. 
(Howlett and Bland, Ref. 110) 
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Fig. 37 Unsteady lower surface pressure distributions for the NACA 64A010A airfoil at 
M - 0.796 and k - 0.1. 

(Howlett and Bland, Ref. 110). 
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a, deg. ’ a, deg. 

(a) Case 1. M = 0.601; (b) Case 2. M = 0.599; (c) Case 3. M = 0.599; ( d ) Case 5. M = 0.7= 


am = 2.89°; ao = 2.41° am = 3.16°; ao = 4.59° am = 4.86°; ao = 2.44° am - 0.02°; rao - 2.5 


Fig. 38 Comparison of unsteady forces versus angle of attack for cases 1, 2. 3. and 5 for the 
NACA 0012 airfoil at k - 0.081 and Re - 4.8 * 106. (Howletl and Bland. Rel. 110) 



ORIGINAL PAGE IS 
OF POOR QUALITY 




(a) NACA 64A010A 
am = 1°. 


(b) MBB-A3 
am = -0.5°. 


ig. 39 Steady pressure distributions. 
(Bland and Edwards, Ref. 114) 
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Fig. 42 Pitch amplitude effect for NACA 64A010A 
airfoil at M = 0.78, k = 0.15. 

(Bland and Edwards, Ref. 114) 


(a) NACA 64A010A .(b) MB3-A3 

Fig. 40 Lifting pressure due to pitch about 
quarter chord at k = 0.15. 

(Bland and Edwards, Ref. 1 1 4) 



(a) Force coefficients. (b) Moment ccefficien 

Fig. 41 Comparison of unsteady aerodynamic 
forces for NACA 64A010A airfoil at 
M = 0.78 and MBB-A3 airfoil at 
M = 0.79. 

(Bland and Edwards, Ref. 114) 
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Fig. 43 Effect of pulse amplitude on lift coefficient 
for NACA 64A010A airfoil. M = 0.796; 
am = -0.21°. 

(Howlett and Bland, Ref. 110) 
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Fig. 45 Sleady pressure distributions for three Fig. 46 Unsteady frequency response functions lor 

airfoil shapes at am = 0°. three airfoils at am - O'*. 

(Batina, 1 Ref. 108) (Batina, Ref. 100) 
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Fig. 47 Steady pressure distributions for three Fig. 48 Thickness effect on unsteady frequency 

airfoil thicknesses at three transonically response functions for three transonically 

scaled Mach numbers at am = 0°. scaled Mach numbers. 

(Batina. Ref. 108 ) (Batina. Ref. 108 ) 
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Fig. 49 Pitching moment coefficient due to flap 
motion for the NACA 64A006 airfoil. 
(Edwards and Thomas. Ref. 1) 
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Fig. 50 Lift-curve for NACA 0012 airfoil 
computed with FL036 code versus 
angle-of-attack as a function of Mach 
number. 

(Salas and Gumbert, Ref. 116) 
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Fig. 51 Anomalous behavior of potential flow 
solutions for NACA 0012 airfoil at 
M = 0.85, am = 0°. 

(Williams et al.. Ref. 118) 
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Fig. 52 Steady pressure distribution for NACA 
0012 airfoil at M = 0.84, am = 0.25’. 
(Fuglsang and Williams, Ref. 119) 
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(b) Unsteady pressures 


Fig. 53 Comparison of measured and calculated 
steady and unsteady pressures for the 
NASA rectangular supercritical wing. 
(Seidel et al.. Ref. 109) 



Fig. 54 Finite-difference grid in the plane of the 
J wing and near- and far-field views. 
(Bennett et al.. Ref. 120) 



Fig. 55 Measured and calculated steady pressures 
for RAE tailplane model at M = 0.80. 
(Bennett et al. Ref. 120) 
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Fig. 56 Measured and calculated unsteady pressure* 
distributions for RAE tailplane model 
oscillating at approximately 70 Hz, 

M = 0.80, k = 0.490. 

(Bennett et al.. Ref. 120). 







a, deg 

Fig. 57 Total lift coefficient for NACA 0012 
wing with AR = 12 at M = 0.82. 
(Gibbons et al.. Ref. 122). 



Fig. 60 Unsteady sectional lift coefficient on 

canard/wing configurations at M = 0.9, 
am = 2 deg and k = 0.3. 

(Batina, ref. 123) 



Fig. 61 Wing-fuselage modeling for TSD potential 
equation. (Batina, ref. 124) 



Fig. 63 Unsteady sectional coefficients due to wing 
first bending for the RAE wing-fuselage at 
M = 0.91, and aw = af = 1°. 



ytdan 3S Experiment ( M = 0.90 ) 

M = 0.90 O Upper surface 

M = 0.91 □ lower surface 



Fig. 62 Comparison between XTRAN3S and experimental wing steady pressure distributions on 
the RAE wing-fuselage at M = 0.90, 0.91, and aw = af = 1°. 

(Batina, ref. 124) 
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64 Comparison of F-5 wing steady pressure Fig. 65 
distributions at M - 0.9 and am = 0°. 



Comparison of F-5 wing upper surface unsteady 
pressure distributions due to wing pitching at M 
am = 0°, ao = 0.109°. and k = 0.137. 




Fig. 67 CAP-TSD modeling of the General Dynamics 
one-ninth scale F-16C aircraft model. 
(Batina, et al., ref. 85) 


ig. 66 Comparison of F-5 wing lower surface unsteady 
pressure distributions due to wing pitching at 
M - 1.1, am = 0°, ao * 0.267°, and k - 0.116. 



Fig. 68 Comparison between CAP-TSD and experimental steady pressure distributions on the 
wing and tail of the F-16C aircraft model at M = 0.9 and am = 2.38°. 

(Batina, et al., ref. 85) 
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Fig. 69 CAP-TSD unsteady pressure distributions on the upper surfaces of the wing and tail of 
the F-16C aircraft model due to complete airplane rigid pitching at M = 0.9, 
am = 2.38°, ao = 0.5°, and k = 0.1. 

(Batina, et al., ref. 85) 



Fig. 70 Planview of 45° sweptback flutter model wing. 
(Cunningham et al., ref. 127) 
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75 Pressure distribution vs. esperiment for 
positive angles of pitching cycle; NACA 0012, 
M = 0.755, k = 0.0814, am = 0.016°, 
ao = 2.51°. 

(Anderson' et al.. Ref. 81) 




. 76 Unsteady forces and moments; NACA 0012, 
Mo. - 0.755, k - 0.0814, am = 0.016°, 
ai = 2.51°. 

(Anderson et al. f Ref. 81) 
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Fig. 77 Lift and pitch-moment coefficients. 

NACA 0012, M~ = 0.755, a 0 = 2.51°, 
k = 0.0814, At = 0.005. 

(Kandil and Chuang, ref. 146) 





Fig. 78 Lift and pitching-moment coefficients, Rg. 
NACA 0012, Moo = 0.6, am = 4.86°, 

ao = 2.44°, k = 0.081, At = 0.01. 

(Kandil and Chuang, ref. 146) 
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79 Summary of surface pressure on a rolling delta 
wing during one cycle of periodic response for 
Moo = 2, a = 10°, © = 0.35, k = 1.337, 9max = If 
(Kandil and Chuang, ref. 146) 








NACA 64A010 

M^-0.8. R«- 12x10® K-0.2 

• UPPER SURFACE 
<* LOWER SURFACE 

FULL 

-- THIN-LAYER 


EXPERIMENT 


I NAVIER STOKES 



t. 


]. 82 Computed variation of aileron angle with time 
showing buzz condition. 

(Steger and Bailey, ref. 25) 
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Fig. 84 Comparison of calculated and measured 
pressure coefficients- for NACA 64A010A 
airfoil at shock-induced separation condition, 
a = 4° + 1° cos <ot. 


(Chyu and Davis, ref. 27) 




Fig. 83 Mean and first harmonic complex components of pressure coefficients: 
a - 0° + 1° cos cot. (Chyu and Davis, ref. 27) 
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a) lift coefficient* 



Fig. 86 Eighteen percent thick circular arc airfc : 
experimental conditions. 

(McDevitt, ref. 4-) 
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S) moment coefficient.' 

Fig. 85 Comparison of unsteady forces versus 

angle of attack for the NACA 0012 airfoil at 
M = 0.599, am = 4.86°, ao = 2.44°, Re * 4.8 x 
(Rumsey and Anderson, ref. 78) 
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Fig. 87 Steady and oscillatory flow domains for 
18 percent thick circular arc airfoil 
(Levy, ref. 45) 


[ j TYPE A □ TYPEC 


END OP UNSTEADY FLOW, 



Fig. 88 Reduced frequency of periodic oscillations on 
an 18 percent thick circular arc airfoil for 
varying angle of attack. 

(McDevitt, ref. 44) 



Fig. 89 Computed Mach contours for 18 percent 
circular arc airfoil, M = 0.72, 0.75. 
and 0.78, Re c = 11 x 106. 

(Levy, ref. 45) 




Fig. 90 Calculated periodic oscillation for 18% biconvex airfoil with implicit thin-layer 
Navier-Stokes code, M « 0.78, Re * 11 x 10 6 , k = 0.406. 

(Edwards and Thomas, ref. 1) 
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Fig. 91 Calculated pitching moment coefficient for 18% biconvex airfoil with non-isentropic 
TSDcode. 

(Edwards and Thomas, ref. 1) 










